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?n*t5eTOurse  of  this  research  we  have  conducted  investigations  of  turbulent  mixing,  chemical  reaction  and 
combustion  processes  in  turbulent,  subsonic  and  supersonic  flows.  The  program  was  comprised  of  several 
parts:  (i)  an  experimental  effort,  (ii)  a  theoretical  and  numerical  simulation  effort,  and  (iii)  a  development 
of  instrumentation  and  diagnostics;  an  upgrade  of  flow  and  combustion  facilities;  and  the  development 
of  data-acquisition  sub-systems.  We  have  carried  out  a  series  of  theoretical  and  experimental  studies  of 
turbulent  mixing  in  primarily  in  two,  well-defined,  fundamentally  important  flow  fields:  free-shear  layers 
and  axisymmetric  jets.  To  elucidate  molecular  transport  effects,  experiments  and  theory  have  dealt  with 
both  reacting  and  non-reacting  flows  of  liquids  and  gases,  in  fully-developed  turbulent  flows.  A  criterion 
for  fully-developed  turbulence  was  recently  developed.  The  computational  studies  are,  at  present,  focused 
at  fundamental  formulation  and  implementation  issues  pertaining  to  the  simulation  of  both  compressible 
and  incompressible  flows  characterized  by  strong  fronts,  such  as  shock  waves  and  flames.  In  our  diagnostic 
development  efforts  we  have  improved  the  signal-to-noise  ratio  of  flow  images,  in  both  ga.s-  and  liquid-phase 
flows,  as  well  as  continued  with  the  development  of  data-acquisition  electronics  to  meet  very  high-speed,  high- 
volume  data  requirements;  the  acquisition  of  single,  or  pairs,  of  two-dimensional  images  in  rapid  succession; 
and  the  acquisition  of  data  from  arrays  of  supersonic  flow  sensors. 
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Abstract 


The  purpose  of  this  research  is  to  conduct  fundamentaJ  investigations  of  tur¬ 
bulent  mixing,  chemical  reaction  and  combustion  processes  in  tiirbulent,  subsonic 
and  supersonic  flows.  The  program  during  this  reporting  period  was  comprised  of 
several  parts; 

a.  201  experimental  effort, 

b.  a  numerical  simulation  effort, 

and 

c.  an  effort  to  develop  instrumentation  and  diagnostics;  flow  and  combustion 
facilities;  and  data-acquisition  systems. 

The  latter  as  dictated  by  the  specific  needs  of  the  experimental  pao-t  of  the  program. 

Our  approach  in  this  research  has  been  to  carry  out  a  series  of  detailed  theoret¬ 
ical  and  experimental  studies  of  turbulent  mixing  in  primarily  in  two,  well-defined, 
fundamentally  important  flow  fields;  free-shear  layers  and  axisymmetric  jets.  To 
elucidate  molecular  transport  effects,  experiments  and  theory  concern  themselves 
with  both  reacting  and  non-reacting  flows  of  liquids  and  gases,  in  fully-developed 
txirbulent  flows,  t.e.,  in  moderate  to  high  Reynolds  number  flows.  A  criterion  for 
fully-developed  turbulence  was  recently  developed  and  will  be  presented  below. 

The  computational  studies  are,  at  present,  focused  at  fimdamental  formula¬ 
tion  and  implementation  issues  pertsuning  to  the  computational  simulation  of  both 
compressible  and  incompressible  flows  characterized  by  strong  fronts,  such  as  shock 
waves  and  flames. 

Our  diagnostic  development  efforts  have  recently  been  focused  on  improving 
the  signal-to-noise  ratio  of  flow  images,  in  both  gas-  Jind  liquid-phase  flows,  as  well 
as  the  continuing  development  of  data-acquisition  electronics  to  meet  very  high¬ 
speed,  high-volume  data  requirements;  the  acquisition  of  single,  or  pairs,  of  two- 
dimensional  images  in  rapid  succession;  and  the  acquisition  of  data  from  arrays  of 
supersonic  flow  sensors. 


1.  Introduction 


Progress  made  under  the  sponsorship  of  this  Grant,  for  the  period  ending  31 
May  1993,  has  taken  place  in  several  focus  areas.  This  is  discussed  in  the  cor¬ 
responding  sections  below.  Additional  documentation  can  be  foimd  in  the  recent 
reports  and  pubhcations  included  as  appendices  to  this  report,  as  well  as  in  the 
reference  list  (Sec.  8).  Additional  copies  of  any  of  this  material  are  available  on 
request. 

Work  in  preliminziry  stages,  at  this  writing,  eind  not  documented  in  the  discus¬ 
sion  below,  includes: 

a.  progress  and  advances  in  signal-to-noise  ratio  in  gas-phase,  supersonic-flow 
imaging;t 

b.  calculations  and  design  for  an  upgrade  of  the  Supersonic  Shear  Layer  (S^L) 
combustion  facility  scheduled  during  the  next  twelve  months; 

c.  an  analytical  effort  to  understand  mixed  (super-/sub-sonic)  flow; 

euid 

d.  advances  in  computer-aided-design  and  fabrication  of  electronic  circuitry'. 


Parts  of  this  effort  were  cosponsored  by  sm  ARPA/ONR  contract,  the  Gels 
Research  Institute,  and  JPL. 


g 


f  Related  work  in  liquid-phase  jets  is  discussed  in  Sec.  4.2. 
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2.  A  criterion  for  fully-developed  turbulence 

Recent  data  on  turbulent  mixing  suggest  that  the  mixing  transition,  previously 
documented  to  occm:  in  shear  layers,  also  occurs  in  the  fax  field  (x/d  >  30)  of 
turbident  jets  jets,  as  well  as  other  flows.  Specifically,  a  transition  to  a  more  well- 
mixed  state  has  been  observed  experimentally  in, 

a.  pipe  flow, 

b.  in  the  near-field  (x/d  <  15)  of  roimd  jets, 

c.  in  cylinder  wakes, 

d.  in  grid  turbulence, 

e.  in  thermal  convection, 

f.  in  Couette- Taylor  flow  (between  concentric  cylinders), 
as  well  as 

g.  in  recent  numerical  simulations  of  turbulence  in  a  spatially  periodic  cube. 

The  resulting,  fully-developed  turbulent  flow  requires  a  minimum  Reynolds  number 
of  Re  «  10^,  or  a  minimum  Taylor  Reynolds  number  of  Rer  «  10^,  to  be  sustained. 
The  transition  itself  must  be  regarded  as  a  more-or-less  universal  phenomenon. 

Turbulent  mixing  in  this  fully-developed  state  does  not  appear  to  be  universal, 
however,  with  a  qualitatively  different  behavior  between  shear  layers  and  jets.  See 
jilso  discussion  in  Miller  &  Dimotakis  (1992).^  These  observations  were  presented 
at  the  Turbulence  Symposium  in  honor  of  W.  C.  Reynolds’s  60‘**  birthday,  22-23 
March  1993  (Monterey,  CA). 

A  GALCIT  report  doc\imenting  this  discussion  is  included  as  part  of  this  report 
as  Appendix  B. 


1  Included  in  this  report  as  Appendix  A. 
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3.  Mixing  and  combustion  in  supersonic,  turbulent  shear  layers 

Several  “flip”  (Koochesfahani  &  Dimotakis  1986),  chemically-reacting,  super¬ 
sonic  sheax-layer  flow  experiments  have  been  performed  in  the  last  year.  These  are 
part  of  an  attempt  to  separate  Reynolds  number  effects  from  compressibility  (Mach 
number)  effects.  In  om  previous  compressible  work  (Hall  ei  al.  1991),  Reynolds 
niunber  effects  were  acknowledged,  but  could  not  be  specifically  addressed. 

The  recent  res\ilts  indicate  that  Reynolds  niunber  effects  are  the  probable  pri- 
meiry  cause  of  lower  mixing  (rate),  at  moderate  compiessibility,  relative  to  that 
documented  on  the  basis  of  previous  data  in  incompressible  reacting  shear  layers, 
at  lower  Reynolds  numbers,  as  opposed  to  compressibility.  This  is  an  important 
issue,  not  only  theoretically  but  also  from  zm  applications  point  of  view.^ 

Previously  available,  chemically-reacting  flow  data  fall  in  the  Reynolds  number 
range  of  10'*  to  10®.  The  minimum  Reynolds  number  for  the  compressible  cases 
is  aroimd  10®.  Important  constraints,  requiring  that  chemical  reactions  proceed  at 
atmospheric  pressiure  to  ensure  that  they  are  in  the  fast-kinetic  regime,  preclude 
compressible  runs  to  be  performed  at  Reynolds  numbers  similar  to  those  of  the 
incompressible  data.  It  is  possible  to  extrapolate  the  incompressible  data  to  these 
higher  Reynolds  numbers,  but  the  uncertainty  of  these  estimates  is  quite  high. 

In  Hall  et  al.  (1991),  the  primary  compressibility  effect  observed  was  between 
the  medium  and  high  compressibility  cases.  Those  experiments,  however,  could  not 
distinguish  between  Reynolds  number  and  Mach  number  effects  and,  as  a  conse¬ 
quence,  it  was  not  possible  to  mslce  definitive  statements  regarding  compressibility 
effects,  at  mediiun  compressibility. 

Figure  1  depicts  the  normalized  product  thicknesses  at  low  stoichiometric  ra¬ 
tios,  <f>,  for  the  compressible  flow  cases  of  Hall  et  al.  (1991)  as  well  as  the  incom¬ 
pressible  flow  ceises  of  Mungal  &  Dimotakis  (1984),  Mungal  et  al.  (1985),  and  FVieler 
(1992).*  The  gap  in  Reynolds  number,  between  the  compressible  and  incompress¬ 
ible  flow  data  is  clear.  The  recent  flip  experiment  pairs  were  designed  to  reduce  the 
uncertainty  in  the  Reynolds  number  effects  at  moderate  compressibility  (Me  «  0.6). 

The  flip  experiment  pairs  were  mewle  with  H2/NO  and  F2  reactants  in  the 


^  D.  Bushnell  (private  communication). 

*  See  also  Dimotakis  (1993),  included  here  as  Appendix  B,  for  additional  discussion. 


Fig.  1  Normalized  product  thickness  for  low  Squares  and  circles:  Mungal  Di- 
motakis  (1984);  vertical  bar:  Mungal  et  al.  (1985);  triangles:  Frieler  (1992); 
arrow  and  diamond:  HaU  et  al  (1991). 

high-  and  low-speed  streams,  respectively,  with  mixtures  of  He  and  Ar  as  diluents. 
Separate  chemical  kinetic  numerical  simulations  and  experimental  investigations 
were  conducted  to  ascertain  that  the  data  were  acquired  in  the  kinetically-feist 
regime.  Fotir  conditions  were  arranged  in  two  pairs.  In  eaich  pwr,  the  compressibility 
(Mach  munber),  density,  velocity  and  sjjecific  heat  ratios  were  kept  fixed.  Within 
each  pair,  the  two  Reynolds  niimbers  were  separated  by  a  factor  of  about  2.5,  t.e., 
nearly  a  half  decade  in  Reynolds  number.  These  experiments  were  designed  such 
that  Reynolds  number  would  be  the  only  cause  of  changes  in  the  chemical  product 
and  mixed  fluid  thicknesses. 

Our  preliminary  recent  results  can  be  compared  with  each  other,  as  well  as 
with  the  existing  body  of  data,  both  compressible  and  incompressible.  In  this  com¬ 
parison,  the  other  differences  in  the  data  must  be  noted.  Most  of  the  incompressible 
flow  pair  data  were  recorded  at  stoichiometric  mixture  ratios  of  ^  =  8,1/8.  The 
flows  at  medium  compressibility  were  nm  with  ^  =  4,1/4.  The  best  assessment 
of  the  difference  can  be  made  in  terms  of  the  two  incompressible  flow  points  from 
Mungal  &  Dimotakis.  These  are  runs  at  identical  conditions  except  that  ^  =  8,1/8 
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for  the  squares,  and  </»  =  4, 1/4  for  the  circles. 

Another  issue  to  consider  is  that  different  data-processing  used  on  the  different 
sets  of  data  might  result  in  differences  between  sets  of  data.  To  this  end,  the  raw 
temperature  measurements  from  Hall  et  al.  and  Frieler  were  reprocessed  using  the 
same  procedures.  Comparisons  of  these  three  sets  of  points  will  have  no  systematic 
error. 

The  normalized  product  thicknesses  for  low  <f>'s  is  shown  in  Fig.  2.  No  evidence 
of  a  Reynolds  number  dependence  is  discernible  within  the  Reynolds  number  span  of 
the  new  (caret  and  star)  data.  This  suggests  that  any  differences  between  compress¬ 
ible  flows  at  these  Reynolds  numbers  is  not  due  to  a  change  in  Reynolds  number. 
There  is  little  difference  between  the  product  thicknesses  extrapolated  on  the  basis 
of  the  incompressible  flow  experiments  and  the  medium  compressibility  cases.  This 
implies  that  the  reduction  of  product  thickness  for  the  compressible  czises,  relative 
to  the  previous,  incompressible  flow,  experiments,  is  due  to  differences  in  Reynolds 
number  and  not  to  compressibility. 


log  (Reg^) 


Fig.  2  Normalized  product  thickness  for  low  <^.  Square  and  circle:  Mimgal  &  Di- 
motakis  (1984);  vertical  bars:  Mungal  et  al.  (1985);  triangles:  Frieler  (1992); 
arrow  and  diamond:  Hall  et  al.  (1991);  carets  and  stars:  recent  experiments. 
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The  normalized  product  thicknesses  derived  from  flows  at  high  <t>  are  shown 
in  Fig.  3.  No  Reynolds  number  effect  is  discernible  in  the  four  new  data  points. 
Differences  between  the  compressible  experiments  axe  probably  not  attributable 
to  the  difference  in  Reynolds  numbers.  The  product  thicknesses  from  the  highest 
Reynolds  munber  incompressible  flow  data  are  comparable  to  the  medimn  com¬ 
pressibility  cases.  Again,  the  implication  is  that  the  effect  of  compressibility  on 
molecular  mixing  is  slight,  if  any,  at  medium  compressibihty. 


4.0  4.5  5.0  5.5  6.0  6.5  7.0 


log  (Reg.^) 

Fig.  3  Normalized  product  thickness  for  high  (j>.  Squeu’e  and  circle:  Mungal  k. 
Dimotakis  (1984);  triangles:  Frieler  (1992);  arrow  and  dijunond:  Hall  et  al. 
(1991);  carets  and  stars:  recent  experiments. 


The  normalized  mixed  fluid  fraction  data,  derived  form  the  high-  and  low-^ 
“flip”  flow  pairs  are  shown  in  Fig.  4.  The  four  new  points  provide  no  clear  evidence 
of  a  Reynolds  number  dependence  for  the  clustered,  compressible  flow  data.  In 
comparisons  with  the  incompressible  data,  it  is  important  to  note  the  two  points 
from  Mung£tl  k  Dimotakis.  As  was  noted  for  the  product  thickness  results,  it  seems 
that  there  is  no  real  reduction  in  mixing  due  to  compressibility. 

Our  recent  preliminary  data  indicate  that  no  strong  dependence  on  Reynolds 
number  exists  between  the  new,  low  juid  high  Reynolds  number,  compressible  flow 
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Fig.  4  Normalized  mixed  fluid  thickness.  Squares  and  circles:  Mimgal  &  Dimotakis 
(1984);  triangles:  Prieler  (1992);  arrow  and  diamond:  Hall  et  al.  (1991); 
carets  and  stars:  recent  experiments. 

cases.  This  indicates  that  differences  between  compressible  flows,  as  with  the  two 
cases  of  Hall  et  o/.,  would  not  be  due  to  Reynolds  number.  Ebctrapolation  of  the 
incompressible  results  produces  a  reasonable  estimate  of  the  mixed  fluid  thickness 
values  of  the  new  data.  This  implies  that  the  documented  reduction  in  the  mixed 
fluid  fraction,  at  moderate  compressibility,  is  due  to  the  higher  Reynolds  numbers 
and  not  the  higher  Mach  niunber  in  the  latter  flows. 

We  appreciate  that  this  conclusic,.'  must  be  regarded  as  preliminary,  in  view 
of  the  data-less  decade,  or  so,  in  Reynolds  munber  across  which  the  extrapolation 
has  been  made.  We  hope  to  be  in  a  position  to  swldress  this  deficiency  in  our  futme 
work. 

This  work  is  part  of  the  Ph.D.  research  effort  of  Chris  Bond. 
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4.  Mixing  and  combustion  in  turbulent  jets 

The  research  effort  on  turbulent  jet  mixing,  involving  both  the  gas-phase,  chem¬ 
ically  reacting  and  liquid-phase,  non-reacting  jet  investigations,  has  been  cospon¬ 
sored  by  the  Gzis  Research  Institute,  GRI  Contract  No.  5087-260-1467. 


4.1  Gas-phase  chemically- reacting  jets 

Previous  experiments  conducted  in  the  High  Pressinre  Reacting  Vessel  (HPRV) 
have  yielded  significant  information  on  the  heat  relezise  distribution  in  a  chemically¬ 
reacting,  axisymmetric  jet.  Those  experiments  docmnented  the  flzune  length  de¬ 
pendence  on  Reynolds  number,  in  the  range  1.0  x  10“*  <  Re  <  15  x  10^,  as  well  as 
the  existence  of  a  Reynolds-number-dependent  mixing  virtual  origin.  In  addition, 
Damkholer  number  effects  on  mixing  and  the  heat  release  distribution  were  inves¬ 
tigated.  See  Gilbrech  (1991)  Dimotakis  et  al.  (1992),  and  Dimotakis  (1993),*  for  a 
more  complete  discussion. 

In  order  to  chsuracterize  flame  length  behavior,  these  investigations  were  con¬ 
ducted  on  a  jet  in  a  purely  momentum-dominated  regime.  Below  are  presented 
the  preliminary  results  of  an  investigation  conducted  on  turbulent  non-premixed 
reacting  jets  m  a  flow  regime  that  cannot  be  regarded  as  in  the  purely  momentum- 
dominated  regime.  New  and  imsuspected  results  were  fovmd  in  both  regions. 

Specifically,  in  the  region  before  the  flame  tip,  there  is  experimental  proof  that 
the  heat  release  occurring  in  the  reacting  zone  of  the  flow  enhances  not  only  the 
entrsiinment,  as  had  been  know  for  some  time,  but  also  the  mixing.  In  the  post¬ 
flame  region,  the  line-integrated  temperatiue  field  is  characterized  by  a  smooth 
dependence  on  the  adiabatic  flame  temperature  rise,  for  a  range  of  adiabatic  flame 
temperatures. 

This  pzirt  of  the  effort  was  principally  undertaken  by  Dr.  Dominique  Fourguette. 


* 


Included  in  this  report  as  Apf>endix  B. 
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4.1.1  Heat  release  effects  in  the  flame  region 

A  detailed  description  of  the  experimental  facility  emd  the  diagnostics  can  be 
foimd  in  the  Gilbrech  (1991)  Ph.D.  thesis,  as  well  as  in  Gilbrech  &  Dimotakis  (1992). 
A  virtually  unconfined,  chemically-reacting  axisymmetric  jet  (2.5  mm  nozzle  diam¬ 
eter)  of  F2,  diluted  in  N2,  is  discharged  into  a  vessel  containing  NO,  also  diluted 
in  N2.  The  adiabatic  flame  temperature  rise,  ATf,  and  the  stoichiometric  mixture 
ratio,  (f>,  axe  set  by  the  relative  concentration  of  the  reactants.  The  nozzle-exit 
Reynolds  number  can  rzinge  over  more  than  an  order  of  magnitude  by  changing 
the  HPRV  ambient  pressure,  thus  not  modifying  other  flow  conditions.  The  ex¬ 
perimental  results  are  in  the  form  of  time-averaged,  line-integrated  temperature 
measurements  both  in  the  flame  region  and  in  the  post-flame  region. 


Fig.  5  Normalized  product  thickness  for  a  purely  and  non-purely  momentum-driven 
flame,  <^  =  14  and  Re  =  1.0  x  10^. 
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Figiire  5  shows  the  line-integrated  temperature  meaisurement  performed  in  the 
purely  momentum-dominated  regime  and  in  the  case  where  heat  release  signifi¬ 
cantly  affects  the  temperature  distribution.  Both  axes  are  normalized  by  the  flzune 
length  measured  in  the  momentum-dominated  regime.  In  the  purely  momentum- 
dominated  czise  (very  low  heat  release),  ATf  =  13.7  K,  the  heat  release  distribution 
is  a  logarithmic  function  of  the  axial  distzmce  up  to  the  fleime  tip  and  is  indepen¬ 
dent  of  the  axial  distance  in  the  post-flame  region  (Gilbrech  1991).  The  logaxithmic 
dependence  is  discussed  in  Dimotakis  (1993),  which  is  appended  to  this  report  (Ap¬ 
pendix  B). 


Fig.  6  Adiabatic  fizune  temperature  scan,  <f>  =  lA  and  Re  =  1.0  x  10^. 


At  higher  adiabatic  flame  temperatxires  (c/.  Fig.  5  data  at  ATf  =  101 K),  the 
line-integrated  measurements  reveal  a  larger  heat  release  in  the  flame  region,  as 
well  as  a  significant  excess  temperature  overshoot  in  the  region  of  the  flame  tip. 
Higher  entrmnment  (of  cold  reservoir  fluid)  aJone  would  result  in  the  same,  or  lower, 
temperatures  in  the  flame  region,  and  therefore  cannot  account  for  the  larger  heat 
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release.  The  higher  line-integrals  of  the  temperature  rise  observed  represent  a  new 
result  and  indicate  an  enhanced  mixing  rate. 


4.1.2  Post-flame  region  temperature  fleld 

The  series  of  rtms  shown  in  Fig.  6  revezd  the  existence  of  a  point  Xc  the 
post-flame  region  where  the  line-integrated  temperature  is  independent  of  ATf .  For 
the  case  <{>  =  14,  this  is  located  at  Xc/Lfo  =  1.96.  In  other  words,  at  a  distance 
downstream  nearly  twice  the  (reference)  flame  length,  Lfo,  i-c.,  the  flame-length 
measured  in  the  (purely)  momentum-dominated  case.  The  relation,  Xc/I»fo  ^  2 
holds  for  all  ^’s  measured,  as  can  be  seen  in  the  following  table.  These  data  also 
suggest  that  the  Re3Tiolds  number  dependence  of  the  normalized  cross-over  p>oint, 
is  carried  by  L{q. 


Xc 

4> 

Re 

XcILfo 

10 

2.0  X  10^ 

2.03 

10 

4.0  X  10-* 

1.94 

12 

1.0  X  10^ 

1.96 

14 

1.0  X  10^ 

1.96 

17.9 

iO  X  10^ 

2.11 

The  consistency  in  the  normalized  cross-over  point,  Xc/Lfo,  provides  further 
credence  to  the  method  utilized  to  estimate  flame  length  and  that  this  estimate  can 
provide  a  useful  scaling  length.  Additional  investigations  are  tmderway  to  probe 
further  into  the  mixing  and  entrainment  mechanisms. 


4.2  Liquid-phase  turbulent  jet  mixing 


In  our  investigation  of  turbulent  mixing  and  scalar  interface  geometry  in  liquid- 
phase  jets,  we  have  obtained  two-dimensional  image  data,  in  a  plane  perpendicular 
to  the  jet  axis,  using  laser-induced  fluorescence  techniques.  High-resolution  (1,024  x 
1,024),  high  signal-to-noise  ratio  (~  250  :  1)  digital  images  over  the  full  transverse 
extent,  in  the  far  field  of  the  jet  interior,  have  been  recorded  with  a  low-noise, 
cryogenically  cooled,  CCD  camera.  In  work  to  date,  such  data  have  been  acquired 
in  the  Reynolds  number  range  3.8  x  10^  <  Re  <  18  x  10^.  For  these  data,  the 
imaging  station  is  at  ar/d  =  300.  Each  of  the  10®  pixels  in  the  CCD-recording  eirray 
is  calibrated  and  normalized  using  ensemble  averages  of  severail  backgroimd  and 
illumination  images.  This  allows  the  instantaneous,  local  jet-fluid  concentration 
field,  c(y,  z,  /  ),  to  be  measmed  at  the  fixed  x-station  imaged  on  the  euray. 

Our  prehminary  data  confirm  that  this  range  of  Reynolds  numbers  spans  a 
mixing  transition  in  the  far  field  of  turbulent  jets,  as  also  documented  elsewhere  in 
this  report  (c/.  Sec.  2  and  Appendix  B  of  this  report).  In  work  in  progress,  we  are 
attempting  to  quantify  this  pre-  and  post-transition  turbulent  mixing  state  using  a 
variety  of  measures. 


5.  Analytical  and  computational  effort 

Our  work  on  the  use  of  Lagrangian  methods  to  compute  unsteady,  1-D  gasdy- 
neimic  flows,  with  and  without  chemical  reactions,  was  recently  pubhshed  (Lappas 
et  al.  1993).  A  reprint  is  included  as  part  of  this  report  in  Appendix  C. 


Fig.  7  Regular  shock  reflection.  Pressure  contour  levels  at  time  t  =  5  (30  contour 
levels  for  0.5  <  p  <  3.5).  Resolution:  60  x  20  cells. 


In  more  recent  work,  a  new  approach  for  computing  multidimensional  flows 
has  been  developed  as  part  of  our  on-going  anaJytical  and  computational  effort. 
The  method  of  characteristics,  which  has  been  successful  for  hyperbohc  systems 
in  two  independent  variables,  has  been  extended  to  the  general  case  of  several 
independent  variables  for  the  particular  case  of  gas  dynamics.  This  new  approach 
was  used  to  develop  a  numerical  scheme  capable  of  computing  multidimensional 
flows  without  the  zirbitrary  dimensional  splitting  that  conventional  methods  employ. 
See  discussion  in  Colella  (1990),  for  example.  A  brief  description  of  the  new  method 
was  documented  in  a  recent  report  (Dimotakis  ei  al.  1992). 


.0  .5  1.0  1.5  2.0  2.5  3.0 

Fig.  8a  Double  Mach  reflection.  Density  contour  levels  at  /  =  0.20  (30  contour 
levels  for  1.7  <  p  <  18.5).  Resolution:  120  x  30  cells. 


1.0 

.8 

.6 


Fig.  8b  Double  Mach  reflection.  Density  contour  levels  at  <  =  0.20  (30  contour 
levels  in  the  range  1.73  <  u  <  21.0).  Resolution:  240  x  60  cells. 


Part  of  our  current  effort  has  focused  on  testing  this  partictdar  method  and 
comparing  it  with  other  conventional  schemes.  We  have  developed  our  own  version 
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of  the  Godunov  scheme  known  as  MUSCL,  which  is  used  for  comparison  purposes, 
and  we  have  zdso  compared  omr  results  with  those  of  other  simvilations  documented 
in  the  literature.  In  pzirticular,  test  cases,  such  as  those  presented  by  Woodward 
and  Colella  (1984),  were  computed  with  very  good  results.  These  cases  include  a 
standzird  oblique  shock  reflection  and  a  more  interesting  case  of  a  Mewrh  reflection. 
See  Figs.  7,  emd  8a,b.  These  results  are  discussed  in  Lappas  (1993). 


.0  .5  1.0  1.5  2.0  2.5  3.0 

Fig.  9a  Inviscid  shear  layer.  Pressure  contour  levels  at  /  =  5.0  (20  contour  levels 
for  0.98  <  p  <  1.52).  Resolution:  200  x  100  cells. 


In  addition  to  validation  nms,  we  are  currently  computing  a  variety  of  flows 
aimed  at  imderstemding  the  physics  of  the  interaction  of  shock  waves  with  a  shear 
layer,  or  contact  discontinuity.  The  computations  are  helping  us  imderstand  the 
often  complicated  wave  patterns  that  arise  in  such  flows.  This  effort  is  helped  by 
the  adaptive  gridding  capability  that  we  have  developed.  A  hiertirchy  of  overlapping 
grids  cein  be  used  to  compute  different  parts  of  the  flow  with  different  resolution. 
This  is  particularly  important  when  computing  regions  of  wave  interaction  where  a 
higher  resolution  may  be  needed.  See  Figs.  9a,b,c. 
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6.  Diagnostics,  Instrumentation,  &  Experimental  Techniques 

Some  of  the  work  performed  as  part  of  this  part  of  the  effort  is  outlined  below. 
The  development  of  the  image  correlation  velocimetry  and  dual-image  CCD  tech¬ 
nology,  in  collaboration  with  JPL,  has  been  cosponsored  by  ARPA/ONR  Contract 
N00014-91-J-1968. 


6.1  Image  correlation  velocimetry 

Oiu:  effort  in  Image  Correlation  Velcomitry  continues  to  refine  the  method, 
introduced  in  last  year’s  report,  for  correlating  two  successive  scalstr  images  for  the 
purpose  of  measuring  imaged  fluid  velocity  field  information.  Specifically,  a  method 
has  been  developed  for  deforming,  or  transforming,  one  image  to  another.  Taylor- 
series  expansions  of  the  Lagrangian  displacement  field  are  used,  in  conjunction  with 
an  integral  form  of  the  equations  of  motion,  to  approximate  this  transformation. 
The  method  locally  correlates  images  for  displacements,  rotations,  deformations, 
and  higher  order  displacement  gradient  fields.  An  integral  form  of  the  equations  of 
motion  is  employed  and,  as  a  consequence,  no  spatial  or  temporal  differentiation  of 
the  image  data  is  required  in  estimating  the  displacement  field.  This  improves  our 
ability  to  handle  data  with  finite  signal-to-noise  ratios.  A  significant  addition  to  our 
method,  over  this  past  year,  is  the  application  of  a  global  minimization  procedure 
to  insure  a  global  consistency  in  the  results.  Successive,  two-dimensional  digitaJ 
CCD  images  of  fluid  motion,  mairked  with  a  passive  scalar,  are  used  to  verify  the 
capabilities  of  the  method.  The  utility  of  the  method  is  also  illustrated  using  a  pair 
of  Voyager  2  images  of  Jupiter. 

A  detailed  account  of  this  effort  to  date  can  be  fovmd  in  Tokumaru  &  Dimotakis 
(1992),  which  is  included  as  part  of  this  report  as  Appendix  D. 


6.2  High-speed  data/image  acquisition  developments 

This  part  of  the  effort  was  principally  imdertaken  by  D.  B.  Lang,  in  collabo¬ 
ration  with  P.  E.  Dimotakis,  as  well  as  J.  Janesick,  T.  Elliot,  and  S.  A.  Collins  at 
JPL. 
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6.2.1  High-speed  A/D  system 

We  currently  have  two  operational,  high-speed  A/D  converter  data  acquisition 
boards  for  high-speed  and  high-resolution  digited  image  as  well  as  scalar  data  2u:qui- 
sition.  The  first  board  has  two  12-bit  15  MHz  A/D  converters  and  can  digitize  two 
analog  inputs  at  15  MHz,  or  a  single  input  at  30  MHz.  The  second  board  has  two 
12-bit  20  MHz  A/D  converters  and  can  digitize  two  inputs  at  20  MHz,  or  one  input 
at  40  MHz.  The  two  bards  can  be  used  in  parallel  for  an  aggregate  data  throughput, 
at  this  time,  of  70  x  10®  data  samples/s. 

Each  board  has  32  MB  of  local  buffer  memory  and  plugs  into  a  standard  VXI 
backplane.  A  VMEbus  computer  running  the  OS-9  real-time  operating  system 
controls  the  A/D  units  through  a  VMEl-to-VXI  bus  converter.  We  intend  to  upgrzwie 
this  to  an  embedded  VXI  computer.  This  will  ehminate  the  need  for  the  VME 
backplane  and  will  zdso  permit  the  display  of  acquired  images  on  the  screen  in  real 
time  for  foctising  and  calibration  purposes. 

Each  A/D  subsystem  is  comprised  of  an  analog  section  eind  a  digital  section. 
The  A/D  converter  analog  section  has  the  following  capabilities: 

a.  two  12-bit  20  MHz  A/D  converters  digitize  two  channels  at  20Msamples/s 
each,  or  one  channel  at  40Msamples/s, 

b.  programmable  gain  of  ±0.4,  ±0.8,  or  ±1.6  V  full  scale, 

c.  programmable  offset  of  ±0.5,  ±1.0,  or  ±2.0  V, 

and, 

d.  Automatic  czdibration  of  A/D  zero  offset  and  full  scjJe  voltage. 

The  A/D  converter  digital  section  has  the  following  capabihties: 

a.  32  MB  of  local  buffer  memory  (70MB/s/board  throughput), 

b.  word-packing  lo^c  saves  memory  and  reduces  data  rate  (reduces  12-bit 
data  rate  from  80MB/s  to  60MB/s), 

and, 

c.  flexible  triggering  and  scanning  options,  i.e. 

i.  trigger  via  front  panel  connectors, 

ii.  trigger  via  CCD  controller  board, 
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iii.  trigger  via  direct  computer  control, 

iv.  32-bit  real-time  coimter  allows  event  time-stamping, 

V.  up  to  8  scan/frame/calibration  patterns, 

vi.  framing  words  mzu'k  beginning/end  of  scan  or  frame, 

vii.  optional  interrupt  at  beginning/end  of  scan  or  frame, 
and, 

viii.  programmable  delay  line  allows  conversion  timing  adjustments. 


These  first-generation  boards  were  designed  and  fabricated  using  computer- 
controlled  wire- wrapping  technology  to  allow  for  small  modifications,  as  dictated  by 
their  first  use.  The  design  can  now  be  frozen  and  we  are  in  the  process  of  redesigning 
them,  relying  on  multi-layer  printed-circuit  board  technology  to  decrease  unit  costs. 
Following  this  development,  four  additional  channels  will  be  fabricated. 


6.2.2  CCD  controller  board 

We  have  recently  completed  the  design  and  fabrication  of  a  very  flexible,  CCD, 
focal  plane  array  controller  bo^u'd  that  can  control  virtually  euiy  standeu'd  or  custom 
CCD  image  array,  including  the  dual-image  Mach-CCD,  presently  imder  develop¬ 
ment  in  collaboration  with  James  Janesick,  Tom  Elliot,  and  Andy  Collins  of  JPL. 
The  board  has  two  microsequencers  that  can  respond  to  external  trigger  events 
and  generate  the  needed  timing  waveforms.  The  pixel  microsequencer  generates 
the  complex  waveforms  needed  by  the  Mach-CCD,  or  any  other  multiphase  clock 
CCD,  and  also  controls  the  A/D  converter  boards.  The  pixel  microsequencer  can 
be  programmed  to  digitize  a  subset  of  the  image  (region  of  interest)  in  order  to 
conserve  memory  and  increase  the  freuning  rate.  The  event  microsequencer  is  used 
to  generate  the  timing  waveforms  for  YAG  lasers,  the  CEunera  shutter,  and  other 
laboratory  devices  that  must  be  synchronized  to  acquire  Both  microsequencers  can 
be  triggered  from  an  external  control  signal,  or  slaved  to  the  other  microsequencer. 
By  way  of  example,  the  real-time  microsequencer  can: 

-  wait  for  an  external  trigger  signal  (e.g.,  the  start  of  a  nm), 

-  fire  the  first  YAG  laser, 

-  signal  the  pixel  microsequencer  to  store  the  first  image. 
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-  fire  the  second  YAG  laser, 
and  finally 

-  signal  the  pixel  microsequencer  to  start  digitizing  the  two  images. 

The  CCD  controller  board  has  the  following  features: 

a.  internal  programmable  clock  of  up  to  40  MHz  (25  ns  timing  resolution), 

b.  two  external  timing/pixel  clock  inputs/outputs, 

c.  three  event  trigger  inputs  (two  coax  and  one  fiber  optic), 

d.  seven  event  trigger  outputs  (four  coax  and  three  fiber  optic), 

and, 

e.  one  37-pin  connector  for  the  camera  head  that  provides: 

i.  the  CCD  clock  input, 

ii.  3  CCD  trigger  inputs, 
and 

iii.  14  CCD  timing  outputs. 

The  fiber  optic  input /outputs  allow  connecting  to  noisy  devices  such  as  the  pulsed 
YAG  lasers. 


6.2.3  CCD  camera  head 

We  have  designed  and  are  currently  fabricating  a  general-purpose  use  CCD 
camera  head  that  czm  be  used  for  the  Mach-CCD.  The  first  version  uses  a  JPL 
camera  head  board  modified  to  enable  acquiring  of  two  images  as  closely  as  1  ^s 
apart.  This  first  version  has  two  anzdog  outputs  providing  up  to  0.5  Mpixels/s  each 
(hmited  by  the  JPL  camera  head  board).  This  allows  the  readout  of  a  pair  of 
1,024  X  1,024  images  in  two  seconds. 

A  second  generation  camera  head  board  is  being  designed  to  provide  4  analog 
outputs  allowing  a  ten-fold  increase  to  5  Mpixels/s  each,  for  an  aggregate  rate  of 
20  Mpixels/s.  The  current  Mach-CCD  has  two  analog  outputs  so  it  will  be  possible 
to  acquire  10  pairs  of  1,024  x  1,024  images  per  second. 
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The  latter  is  a  very  significant  development,  and  will  permit  a  much  better 
utilization  of  the  experimental  facilities  used  as  part  of  this  research  effort.  By 
way  of  example,  the  cost  of  some  high  Mach  number,  chemically-reacting  runs  in 
the  supersonic  shear  layer  facihty  can  exceed  $2,000,  each.  At  present,  we  can 
acquire  only  a  single,  high-quality  image  in  such  a  nm.  The  new  gereration  of 
image-acquisition  technology  xmder  development  should  permit  higher  signal-to- 
noise  ratios,  along  with  a  20-fold  increase  in  the  image  framing  rate. 
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M.  D.  Slessor:  Graduate  Research  Assistant,  Aeronautics; 

P.  Svitek:  Member  of  the  Technical  Staff,  Div.  of  Engineering  &  Ap¬ 
plied  Science; 
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Abstract 

Single-point,  jet-fluid  concentration  measiurements  obtained  from  high  Schmidt 
number  (Sc  ~  1.9  x  10’)  turbulent  jets  permit  an  investigation  of  temporal  scalar 
power  spectra,  for  jet  Reynolds  numbers  in  the  range  of  1.25  <  Re  x  10”^  <  7.2 . 
At  intermediate  scales,  we  find  a  spectrum  with  a  logarithmic  derivative  (sloi>e) 
that  is  increasing  with  Reynolds  number,  in  absolute  value,  but  less  than  */3  S’t  the 
highest  Reynolds  number  in  our  experiments.  At  the  smallest  scales,  our  spectra 
exhibit  no  power-law  behavior,  possessing  a  log-normal  region  over  a  range  of 
scales  exceeding  a  factor  of  40,  in  some  cases. 
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1.  Introduction 

We  report  here  on  an  experimental  investigation  of  temporal  scalar  power  spec¬ 
tra  of  round,  high  Schmidt  number,  turbulent  jets.  It  is  part  of  a  larger  effort  to 
better  understand  mixing  in  turbulent  free  shear  flows,  including  jets  and  shear  lay¬ 
ers.  In  these  experiments,  we  examined  the  jet  fluid  concentration  (scalar)  power 
spectra  for  several  reasons.  Spectra  are  sensitive  diagnostics  of  the  flow,  providing 
information  over  a  wide  range  of  scales.  Historically,  they  have  been  the  object  of 
a  great  deal  of  attention,  partially  because  it  is  possible  to  extract  predictions  for 
spectral  slopes  from  various  turbulence  theories  and  models. 

Key  among  these  turbulence  theories  are  the  1941  paper  by  Kolmogorov,^ 
with  implications  for  the  scalar  field  in  the  inertial  rainge  discussed  by  Corrsin,^ 
Oboukhov,^  and,  for  higher  wavenumbers,  the  theory  by  Batchelor.^  See,  for  exam¬ 
ple,  discussions  in  Monin  and  Yaglom,®  as  well  as  in  the  recent  review  by  Gibson.® 
Both  the  Corrsin  and  the  Oboukhov  theories  yield  predictions  of  power-law  spectra 
and  of  the  spectral  power- law  logarithmic  derivative,  or  slope,  as  it  will  be  subse¬ 
quently  referred  to  in  this  paper.  Specifically,  the  Corrsin  and  Oboukhov  theories 
predict  a  scalar  spectrum  proportional  to  in  the  inertial  range,  as  did  the 

1941  Kolmogorov  theory  for  the  energy  spectrum. 

For  energy  spectra,  this  has  been  observed  experimentally  under  many  condi¬ 
tions  (c/.  compilation  of  data  by  Chapman^).  The  situation  is  less  clear  concerning 
scalar  spectra,  with  departures  from  the  predicted  behavior  continuing  to  fuel  de¬ 
bate  about  details  and  refinements  of  the  theory. 

Batchelor^  and  Batchelor  ei  al.^  recognized  that  for  Schmidt,  or  Prandtl,  num¬ 
bers  away  from  unity  there  exists  an  additional,  scalar-diffusion,  scale,  now  com¬ 
monly  referred  to  as  the  Batchelor  scale,  which  admits  a  change  in  the  scalar  sp>ec- 
tral  behavior.  The  Batchelor  theory^  predicted  that  the  scalar  power  spectrum  at 
high  Schmidt  numbers  would  display  &  k~^  dependence  beyond  the  Kolmogorov 
wavenumber,  t.c.,  a  spectral  slope  of  —1.  Measurements  in  the  laboratory  {e.g., 
Gibson  and  Schwarz®)  and  the  ocean  {e.g..  Grant  et  al.^°)  were  subsequently  re¬ 
ported  to  be  in  accord  with  this  prediction. 
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On  the  other  hand,  more  recent  measurements  by  Gargett^^  in  the  ocean  were 
found  not  to  exhibit  a  k~^  spectral  range  (see,  however,  discussion  by  Gibson  in 
Refs.  12  and  6).  The  same  result  was  noted  in  pztssive  scalzu:  mixing  measurements  in 
shear  layers^^  zind  in  sczdar  measurements  in  turbulent  jets.^"*’^®  Specifically,  despite 
adequate  resolution  in  those  experiments,  no  k~^  range  was  foimd  at  high  spatial 
wavenumbers,  or,  to  be  exact,  temporal  frequencies,  adding  to  the  questions  about 
the  imiversality,  if  not  the  validity,^®  of  the  classical  predictions  at  high  Schmidt 
number. 

The  issue  of  spatial,  vs.  temporal,  spectra  should  be  recognized  here.  The 
classical  theories  cited  deal  with  spatial  spectra.  One  could  argue,  therefore,  that 
comparisons  of  measurements  of  temporal  spectra  with  predictions  of  spatial  spectra 
cannot  be  made  directly.  Two  points,  however,  should  be  noted  in  response.  First, 
the  overwhelming  majority  of  experimentally  obtained  spectra  reported  in  accord 
with  the  theoretical  predictions  have,  in  fact,  been  temporal.  Second,  at  least  in 
the  case  of  developing  flows  that  are  not  (statistically)  spatially  homogeneous  over 
the  range  of  spectral  scales  of  interest,  the  notion  of  a  spatial  spectrum  and  the 
assumption  of  a  statistically  spatially  homogeneous  turbulent  field  is  questionable. 
Temporal  spectra,  derived  from  point  measurements,  do  not  have  to  contend  with 
this  issue. 


2.  Experiment 

The  experiments  investigated  the  scalar  (concentration)  field  of  round,  axisym- 
metric,  momentxim-dominated,  txirbulent  jets  issuing  from  contoured  nozzles  into  a 
large,  quiescent  discharge  tank.  The  measurements  were  performed  in  the  far  field, 
on  the  centerline  of  the  jet.  Details  of  the  experimental  apparatus  have  appeared 
previously,^^"^^  so  only  a  brief  overview  will  be  presented  here. 

The  experimental  facility  consists  of  three  major  parts:  the  jet  plenum,  noz¬ 
zle,  and  delivery  system;  a  large  reservoir  that  acts  as  the  discharge  tank;  and 
the  diagnostics,  consisting  of  an  argon-ion  laser,  focusing  optics,  collection  optics, 
detector,  signal-processing  electronics,  and  the  subsequent  data-processing.  The 
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working  flmd  is  water,  and  the  scalar  is  a  laser  dye  (disodivun  fluorescein)  which  is 
homogeneously  premixed  with  the  jet  plenum  fluid.  The  resulting  Schmidt  number, 


(1) 


with  V  the  kinematic  viscosity  of  water  and  ~  5.2  x  10  ®  cm^/s  the  estimated 
aqueous  species  diffiisivity  of  the  fluorescein  dye  (Ref.  18,  p.  280),  is  1.9  x  10^. 


The  jet  flow  was  established  and  mcuntained  by  pressurizing  a  downward  ori¬ 
ented  jet  plenum  with  gas.  Both  sonically-metered  and  blow-down,  nearly  constant 
pressure,  gas  delivery  configurations  were  used.  The  internal  exit  diameter  of  the  jet 
nozzle  is  2.5mm  (0.1  in.).  The  rectangular  discharge  tank  is  square  in  cross-section, 
approximately  2  m  high  and  1  m  on  edge.  The  tank  bottom  is  over  600  nozzle  di¬ 
ameters  downstream.  Large  glass  windows  on  all  four  sides  provided  optical  access 
(see  Ref.  15  for  details). 

The  illumination  source  was  an  argon-ion  laser  (Coherent  Innova  90).  The 
particular  unit  was  custom  selected  for  its  low  AM  noise  figure  — 95dB)  over 
the  frequency  range  of  interest  in  these  experiments.  It  was  operated  at  a  power 
of  3.5  W  in  the  light-regulation  mode.  The  beam  was  spatially  filtered,  expanded, 
collimated,  and  subsequently  focused  to  a  small  waist  located  on  the  centerline  of 
the  jet.  A  low  dye  concentration  was  used  in  the  jet  plenum  (~  10“*M),  with 
correspondingly  substantially  lower  concentrations  at  the  measuring  stations.  A 
more  detailed  discussion  of  this  and  related  issues  may  be  found  in  Refs.  14  and  15. 
The  emitted  fluorescence  intensity  was  then  proportional  to  the  local  scalar  (dye) 
concentration  c(x,  t),  that  was,  in  t\im,  averaged  over  the  extent  of  the  measurement 
volume. 

The  fluorescence  emitted  from  the  measurement  volume  was  collected  through 
a  narrow  slit  spatial  filter.  The  beam  profile  and  the  slit  width  defined  a  small, 
spatially-averaging  volume,  roughly  spherical  in  shape  and  of  extent  (diameter) 
^4  fn  50 /4m,  as  estimated  by  direct  observation  using  a  cathetometer.  We  will 
return  to  this  quantity  later. 
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A  photomultiplier  tube  (RCA  8645)  was  used  to  detect  the  fluorescence  emitted 
from  this  volume.  Its  output  signal  was  amplifled  by  a  custom-designed  low-noise 
transimpedance  amplifier,  low-pass  filtered  using  a  third-order  Butterworth  filter, 
digitized,  sampled  with  some  margin  with  respect  to  the  Nyquist  frequency,  and 
stored  for  subsequent  processing. 

The  measurements  to  be  discussed  here  were  made  in  the  fax  field,  on  the  axis 
of  the  jet,  for  jet  Reynolds  mmibers  in  the  range  of 

1.25  X  10^  <  Re  =  —  <  7.2  x  10^  ,  (2) 

V 

where  uo  is  the  nozzle  exit  velocity,  d  the  nozzle  exit  diameter,  and  u  is  the  kinematic 
viscosity.  Data  were  also  recorded  at  both  lower  and  liigher  Reynolds  numbers.  The 
lower  Reynolds  number  jets,  however,  behaved  substantially  differently,  by  any  of  a 
number  of  criteria,  and  were  not  accepted  as  representing  bona  fide  turbulence.*^  In 
the  other  limit,  the  jet  at  the  higher  Reynolds  number  (Re  =  10.2  x  10^ )  produced  a 
distinct  hissing  sotmd.  This  was  probably  generated  by  the  transient  dilatation  and 
subsequent  oscillations  of  small  air  bubbles  caused  by  the  rapid  reduction  in  pressure 
in  such  bubbles  as  they  exited  the  nozzle,  or  by  cavitation  in  the  jet  near-field  region, 
or  both  (note  that  the  plenum  gauge  pressure  is  quadratic  with  Reynolds  number). 
See  discussion  and  references  in  Ref.  19,  pp.  452-453,  and  Ref.  20,  pp.  205-207,  for 
example.  As  a  result,  the  Re  =  10.2  x  10^  jet  was  exposed  to  different  near-field 
conditions  and  will  not  be  included  in  the  discussion  below.  See  Ref.  15  for  a  further 
documentation  of  the  data. 

Finally,  constraints  dictated  by  resolution  and  statistical  convergence,  vis-d-vis 
total  number  of  large  scale  structures  captiued  and  length-of-nm  considerations, 
led  to  measurement  stations  in  the  range, 

100  <  3  <  305  ,  (3) 

d 


where  x  is  the  distance  from  the  nozzle  exit. 
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3.  Scalar  power  spectrum  estimation 

The  fluorescence  signzd  <f>(t),  representing  the  photon  flux  incident  on  the  pho¬ 
todetector,  is  a  linear  function  of  the  spatied  average  of  the  convected  local  jet-fluid 
concentration  c(x,  t)  over  the  measurement  volume.  It  produces  a  signal  that  can 
be  approximated  by  a  convolution  over  c(t)  =  c(xo,<),  the  jet  fluid  concentration 
at,  say,  the  center  of  the  measurement  volume,  Xq,  i.e., 

<l>(t)  ~  f  h^(t -i*)c(t')dt'  =  hf,(i)®c(t)  .  (4) 

•/  — OO 

In  this  expression,  h,^(i)  models  the  impulse  response  of  the  spatio-temporal  aver¬ 
aging  process,  i.e.,  the  temporal  signal  that  would  be  measured  if  a  spatial  delta 
function  of  dye  was  convected  through  the  measurement  volume  at  the  local  flow 
velocity.  See  Fig.  1. 


C(X,  t) 


Fig.  1  Sketch  of  jet  fluid  concentration  field  c(x,i),  convected  through  the  mea¬ 
surement  volume  of  extent  £o- 

The  fluorescence  output  ^(t),  along  with  fluctuations  contributed  by  the  small 
laser  intensity  fluctuations,  convected  residual  non-uniformities  in  the  jet  plenum 
dye  concentration,  photon  shot  noise,  electronic  noise  generated  by  the  signal¬ 
processing  chain,  etc.,  was  processed  by  the  Butterworth  low-pass  filter  to  produce 
the  total  signal 


s(t)  =  hLp(<)®^(t)  +  n(t)  =  h(t)®c(t)  -f  n(t)  , 


(5) 
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that  was  digitized  and  stored.  In  this  expression,  h{t)  =  /iLp(t)  0  is  the  total 
system  transfer  function  2ind  n(t)  is  the  total,  low-pass-filtered  noise. 

Assuming  that  the  total  system  noise  n(/)  can  be  modeled  as  imcorrelated 
with  the  local  dye  concentration  time  history  c(<),  the  spectrum  S,{u})  of  s{t)  can 
be  expressed  in  terms  of  the  spectrum  5^(u;)  of  <j>{i)  and  the  spectrum  5n(u;)  of  the 
(low-pass  filtered)  noise  n(t),  i.e., 

Ssiu)  ~  S^iu)  -h  S„(u;)  ,  (6) 

where,  from  Eq.  4, 

S,(u,)  5.(0.)  ,  (7) 

with  Ha,{u)  =  ^T{  },  the  Fourier  transform  of  /ia(<).  This  allows  us  to  relate 

the  total  signal  spectrum,  5*(a;),  to  the  desired  scalar  fluctuation  spectrum,  Sc(u;), 
of  c{t),  i.e. 

S,(u)  -  mu)\^Sc(u)  -f  SnH  -  +  S„(i^)  ,  (8) 

where  h(i)  }  =  ffLPCw)  Ha(u;). 

For  these  experiments,  the  knee  of  the  Butterworth  low-pass  filter  was  set  sub¬ 
stantially  higher  than  the  range  of  frequencies  contained  in  5^(u;).  Its  main  purpose 
was  to  bandlimit  the  noise  and  de-alias  the  digitized  measurements,  allowing  the 
noise-floor  to  be  determined,  as  will  be  illustrated  in  the  spectra  presented  below. 
This  is  the  reason  the  modulus  squared  of  /fLp(‘*')j  the  transfer  function  of  the 
low-pass  filter,  can  be  ignored  in  Eqs.  7  and  8,  and  wherever  it  multiplies  5^(u;) 
and  5c(u?). 

Figure  2  illustrates  these  relations  by  comparing  the  spectrum  5«(u;)  of  the 
total  signal  s(t),  i.e.,  fluorescence  ^(t)  plus  noise  n(t),  with  5^(u>),  the  spectrum  of 
the  fluorescence  signal  alone.  The  latter  was  calculated  by  subtracting  the  estimated 
noise  spectrum,  S„(uj),  from  the  total  spectrum  5,(u;).  Recalling  Eq.6,  we  have 


5^(u;)  ~  S,(uj)  -  5n(u;)  . 


(9) 
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The  noise  spectrum  was  assiuned  to  be  white,  as  was  found  to  be  the  case  in  sepa¬ 
rate  measurements  of  this  quantity  (see  also  Ref.  21  for  examples).  Nevertheless,  the 
result  is  not  sensitive  to  the  assumed  shape  of  the  noise  spectrum  at  low  frequencies 
where  S^(u))  dominates.  The  data  processed  to  produce  the  spectra  in  Fig.  2  were 
recorded  at  x/d  =  100,  for  Re  =  1.25  x  10^.  Note  the  high  d3meunic  reinge  of  the 
total  signal  spectrum,  t.c.,  the  (logarithmic)  difference  of  the  low-frequency  power 
to  noise-floor  power.  Note  also  that  the  span  to  one-half  the  (scaled)  sampling  fre¬ 
quency  is  well  beyond  the  noise-cross-over  frequency.  As  can  be  seen,  the  frequency 
extent  of  the  noise  floor  was  substantial. 


lOfliolf  Tj) 


Fig.  2  Sample  spectrum  of  the  total  signal  (solid  line:  fluorescence  -|-  noise)  and 
estimated  fluorescence  spectrum  (dashed  line:  fluorescence),  derived  from 
measurements  at  x/d  =  100,  He  =  1.25  x  10^.  Frequency  scaled  by  Tf(x), 
the  local  large-scale-passage  time. 


The  spectra  in  Fig.  2,  and  throughout  this  paper,  are  normzdized  by  c^,  the 
square  of  the  local  mean  value  of  c(f),  multiplied  by  the  local  large-scale-passage 
time,  Ts(x),  and  plotted  in  terms  of  the  circular  frequency,  /,  scaled  by  rf(x).  In 
these  coordinates,  their  integral  produces  the  normalized  variance,  e.g., 


o'  Jo  I  c‘  Tf  i 


(10) 
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The  local  large-scale-passage  time,  rj,  is  given  by, 


Ts{x) 


^(x) 

«cl(x) 


(11a) 


where 

6{x)  0.41  X  ,  (11b) 


is  the  local  outer  scale  of  the  flow,  here  identifled  with  the  (measured)  mean  trans¬ 
verse  extent  (visual  width)  of  the  conical  region  enveloping  the  jet-fluid  (Ref.  15, 
Appendix  D),  and  Uci(x)  is  the  mean  centerline  velocity.  The  latter  was  estimated 
from  the  relation 


Uci(x) 


Ui 


=  6.2 


X  —  X; 


(lie) 


where  uj  is  the  jet  velocity  and  xj  the  jet  (virtual)  origin,  recommended  by  Chen 
and  Rodi.^^  This  spectrum  and  frequency  scaling  was  foimd  to  produce  similarity 
with  respect  to  the  downstream  coordinate,  x/d,  in  the  analysis  of  scalar  spectra 
meastired  in  gas-phase  jets.^* 


The  spectrum  Sc{u)  of  the  scalar  fluctuations  c{t)  can,  at  least  formally,  be 
estimated  by  solving  Eq.  7,  ».e., 

or  Elq.  8,  yielding  a  result  in  terms  of  experimentally  estimated  quantities,  t.e.. 


5c(w) 


S,{u})  -  Snji^) 


(12b) 


We  should  note,  at  this  point,  that  Eq.  4,  assuming  a  fixed  h^{t),  is  not  a  proper 
equation  for  two  reasons.  First,  the  fluid  velocity  convecting  the  c(x,  t)-field  through 
the  measurement  volmne  is  not  a  constant.  Second,  different  scalar  field  lagrangian 
trajectories  through  the  measmement  volume  sample  chords  of  diffo-ent  sizes  (and 
transit  times)  through  it.  These  two  effects  could,  in  principle,  be  expressed  as  two 
additional  convolutions  over  the  local  velocity  field  and  scalar  paths  through  the 
measurement  volume.  See  Fig.  1. 
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For  frequencies,  however,  corresponding  to  spatial  scales  of  the  order  of,  or  less 
than,  the  extent  of  the  measurement  volume,  where  the  transfer  fimction 

will  have  an  effect  on  the  measurement  of  the  scalar  spectrum,  fractional  fluctua¬ 
tions  in  the  convecting  velocity  are  considerably  smaller  than  fractional  fluctuations 
in  c(x,t),  at  leeist  for  this  high-Schmidt-number  fluid.  They  are  also,  largely,  uncor¬ 
related  with  them.  Fluid  paths  through  the  measurement  volume  are  also  imcorre- 
lated  with  the  passively-convected  scalar  field.  As  a  consequence,  in  estimating  the 
spectrum  in  the  frequency  neighborhood  of 

>  (13) 

and  above,  these  two  effects  do  not  contribute. 


Alternatively,  the  most  generzil  expression  of  the  linear  dependence  of  the  flu¬ 
orescence  signal,  on  the  local  jet-fluid  concentration,  c(t),  is  given  by  Eq.  7. 
While  Eq.  7  follows  from  Eq.  4,  the  converse  is  not  true.  Equation  4  is  more  re¬ 
strictive,  also  requiring  a  definite  phase  relation  between  c(<)  and  ^(<).  Fortunately, 
these  (unknown)  phase  relations  do  not  enter  in  the  relation  between  the  corre¬ 
sponding  spectra.  We  may  conclude  that  Eqs.  7  and  12  represent  valid  relations  for 
the  corresponding  spectra,  even  as  Eq.  4  cannot  be  accepted  as  a  correct  description 
of  the  time-history  of  the  fluorescence  signal 


The  estimation  of  the  scalar  spectrum  5c(u;)  in  the  frequency  rzmge  influenced 
by  the  spatio-temporal  averaging  of  the  measurement  process  also  requires  knowl¬ 
edge  of  the  transfer  function.  This  can  be  estimated,  in  turn,  by  noting  that 

it  is  dominated  by  a  pole  corresponding  to  the  transit  time  of  the  flow  through  the 
measurement  volume,  t.e.. 


with 


1 

1  4-  iu;r»  ’ 


(14a) 


(14b) 


Performing  two  different  experiments,  under  as  identical  flow  conditions  as  was 
feasible,  at  two  different  spatial  resolutions,  we  were  able  to  compare  the  spectra, 
for  p  =  1,2, 


5,,(u.)  =!  |ff,(u,)|' S.(w)  +  5„,(«)  , 


(ISa) 
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corresponding  to  two  different  locations  of  the  dominant  pole,  at,  say,  =  tj  and 
r,  =  T2  >  Ti ,  *.e.,  for  p  =  1,2, 


1 

1  +  tWTp 


(15b) 


If  the  sczdar  spectrum,  ScCw),  could  be  assiuned  to  be  identical  in  the  two  experi¬ 
ments  (the  corresponding  noise  floors  were  determined  separately  in  each  case),  we 
see  that  the  ratio  of  the  two  fluorescence  spectra. 


£^1  ^  S»A^)-SnM 


G{u)\Ti,T2)  , 


would  be  given  by 


Giu}-,Ti,T2)  ~ 


1  -|-  (u»’T2)^ 
1  +  (wTi)^ 


(  1 ,  for  w  <  l/r2  ; 

\  (Ti/Tif  ,  for  u;  >  1/ri  , 


(16a) 


(16b) 


independently  of  the,  as  yet  tmknown,  scalar  spectrum,  5c(u;). 


logjot^T,) 


Fig.  3  Dotted  lines:  Fluorescence  spectra  estimated  from  measurements  at  x/d^ 
100  and  Re  =  1.25  x  10^  at  two  spatial  resolutions.  Solid  line:  Estimated 
concentration  spectrum  (Eq.  12)  at  z/d  =  100  and  He  =  1.25  x  10^.  Dashed 
line:  Estimated  concentration  spectrum  at  x/d  =  305  and  He  =  1.2  x  10^. 
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Two  fluorescence  spectra,  =  5*, (o')  —  S„^{u;)  and  —  Sa,  (w)  -  5„,(a;), 
from  one  such  pair  of  experiments,  at  x/d  =  100  and  J?e  =  1.25  x  10^,  are  plotted 
in  Fig.  3  (dotted  lines).  The  fluorescence  spectrum  with  the  larger  high-frequency 
content  is  the  one  plotted  in  Fig.  2.  Power  spectra  were  computed  numerically 
using  a  power  spectral  density  estimation  methodology  that  hcis  evolved  over  the 
past  ten  yeeirs,  or  so.  A  documentation  of  some  of  its  earlier  features  can  be  foimd 
in  Ref.  24.  In  processing  the  data  in  the  experiments  reported  here,  the  power 
spectral  density  estimation  program  computes  spectra  of  data  files  by  means  of 
FFT  methods,  eind  incorporates  Hanning  windowing,  contiguous  record  overlapping, 
Eind  parabolic  detrending,  among  other  featm-es.  Records  up  to  2^^  points  cem  be 
accommodated.  For  spectra  known  to  be  smooth,  the  program  can  provide  third- 
octave  (~  1/10  decade)  gaussian  filtering,  sampled  at  20  points  per  decade,  to 
produce  the  final  spectra.  This  feature  was  used  for  all  the  spectra  plotted  in 
this  paper,  not  so  much  for  smoothing,  but  to  reduce  the  number  of  points  to  a 
manageable  level  for  plotting  purposes  (note  that  2^®  =  65,536). 


looioft  Tj) 

Fig.  4  Computed  ratio,  G(u;;ti,T2),  of  fluorescence  spectra  at  x/d  =  100  (Eq.  16). 
Circles:  Re  =  1.25x10^.  Squares:  Re  =  2.55x10^.  Solid  line:  Least  squares 
fit  for  Ti  zmd  T2. 
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The  ratios  G(a;;Ti,T2)  for  a  pair  of  spectra,  measured  at  x/d  =  100,  at  Re  = 
1.25  X  10^  (circles)  auid  a  pair  at  Re  =  2.55  x  10^  (squares),  respectively,  are  plotted 
in  Fig.  4.  As  c^ln  be  seen,  the  ratios  of  spectra  measured  at  different  Reynolds 
numbers  are  very  nearly  the  same,  in  accord  with  the  analysis  outlined  above,  even 
though,  as  we  will  see  below,  the  spectra  themselves  are  Re-dependent. 

The  curve  (solid  line)  in  Fig.  4  is  a  least-squares  fit  for  rj  and  r2  to  the  lower 
Reynolds  number  data,  that  were  characterized  by  the  higher  signal-to-noise  ratio, 
in  the  frequency  range  1.0  <  frs  <  2.8.  The  lower  limit  of  the  fit  range  is  chosen 
so  as  to  exclude  (the  small)  run-to-run  variations  at  frequencies  well  below  those 
affected  by  the  spatial  averaging.  The  upper  frequency  limit  is  dictated  by  the  less 
than  unity  signal-to-noise  ratio  at  higher  frequencies  yet  (c/.  Fig.  2).  The  values  for 
Ti  and  T2  estimated  by  this  procedure  were 

2irTi/Ts  4.2  X  10“^  and  2irT2/Ti  ~  6.6  x  10“^  ,  (17) 

respectively.  This  corresponds  to  an  effective  spatial  extent  of  the  measurement 
volume  of  £»  —  69  fim  for  the  smaller  of  the  two  (Eq.  13),  in  reasonable  accord 
with  the  visually  estimated  value  of  ~  50  fim  using  the  cathetometer.  This  value 
was  used  to  calculate  the  concentration  spectrum,  5c(u?),  at  x/d  =  100  and  Re  = 
1.25  X  10^,  plotted  as  the  solid  line  in  Fig.  3.  It  was  computed  from  the  fluorescence 
spectrum  recorded  at  the  higher  resolution,  using  Eq.  12  with  the  estimated  single¬ 
pole  transfer  function  Ra(u;)  of  Eq.  14,  at  =  tj. 

The  effective  pole  locations  for  data  recorded  at  x/d  ■=  305  were  more  diflicult 
to  estimate.  At  x/d  =  305,  the  higher  relative  spatial  resolution  pushed  the  poles 
closer  to  the  noise  cros.s-over  point.  On  the  other  hand,  at  x/d  =  305,  the  (logarith¬ 
mic)  difference  between  the  fluorescence  and  estimated  concentration  spectra  was 
much  smaller  over  the  frequency  range  of  interest.  Figure  5  plots  the  fluorescence 
spectrum  (dotted  line)  at  x/d  =  305  and  Re  =  1.2  x  10^  as  weU  as  the  estimated 
concentration  spectrum  (dashed  line).  As  can  be  seen,  the  effects  of  compensation, 
in  this  case,  are  much  smaller  (c/.  difference  at,  say,  /t«  ^  2.7  in  Figs.  3  and  5). 
The  estimated  concentration  spectrum  at  x/d  =  305  in  Fig.  5  is  the  one  plotted  as 
a  dashed  line  in  Fig.  3. 
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Fig.  5  Dotted  line:  Fluorescence  spectnun.  Dashed  line:  Estimated  concentration 
spectrum.  Data  recorded  at  x/d  =  305  and  Re  =  1.2  x  10^. 

The  jet  fluid  concentration  spectra  to  be  discussed  below  were  all  estimated 
in  this  fashion.  The  values  of  the  transfer  flmction  time  constant  used  in  the 
compensation  calculations,  were  fixed  for  all  the  data  measured  at  each  x/d  axial 
location  (Eq.  17,  for  measmrements  at  x/d  =  100,  depending  on  which  slit  width 
was  used  to  record  the  fluorescence  data).  A  fixed  pair  of  %alues  was  also  used  for 
all  the  data  measured  at  x/d  =  305. 


4.  Results  and  discussion 

The  conspicuous  agreement  between  the  concentration  spectra  at  x/d  —  100 
(solid  line)  and  x/d  =  305  (dashed  line)  in  Fig.  3,  up  to  firequendes  limited  by 
signal-to-noise  ratio  considf  ations,  should  be  noted.  A  similar  independence  of  the 
scaled  spectra  with  dowTisti  oca  location  was  also  fotmd  to  hold  in  gas-phase  jets,^^ 
where  the  relatively  lai^v  uffusion  scales,  at  Sc  »  1,  and  comparable  Reynolds 
numbers  made  it  possible  to  estimate  the  concentration  spectra  with  emmgh  spatial 
resolution  directly,  obviating  the  need  for  the  compensation  scheme  employed  here. 
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In  the  case  of  the  present,  liquid-phase  data,  we  should  entertain  the  possibihty 
that  the  collapse  of  the  spectra  may  be  attributable  to  a  fortuitous  choice  of  the 
compensation  pole.  Some  observations,  however,  are  relevant  here.  First,  the  close 
agreement  between  the  two  spectra  also  holds  in  the  /r^  <  2  frequency  range, 
where  the  effects  of  compensation,  even  for  the  x/d  =  100  data,  are  negligible. 
As  we  will  see  below,  this  will  also  be  found  to  be  the  case  at  higher  Reynolds 
numbers.  Second,  the  value  of  the  pole  represents  a  single  parameter.  In  contrast, 
the  collapse  of  the  spectra  at  the  two  axial  locations,  with  substantially  different 
degrees  of  compensation  in  each  case,  is  over  5  orders  of  magnitude  in  the  power 
spectra.  Conversely,  this  collapse  supports  the  validity  of  the  single-pole  model  for 
the  leading  order  behavior  of  the  measurement  transfer  fimction. 


Fig.  6  Frequency-scaled  concentration  spectra.  Solid  line:  x/d  =  100  and  Re  = 
1.25  X  10^.  Dashed  line:  x/d  =  305  and  Re  =  1.2  x  10^.  Dotted  line: 
reference  line  at  a  Va  slope,  corresponding  to  a  k~^  spectrum. 


It  is  useful  to  plot  the  product  of  the  concentration  spectra  with  as  is 

commonly  done.  A  spectnun  described  by  a  -Va  power-law  yields  a  horizontal 
line  over  the  -5/a  frequency  range  when  plotted  in  this  fashion.  The  product  of 
the  concentration  spectra  with  (/r<)*/’,  derived  from  the  data  at  x/d  =  100  and 
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x/d  =  305,  at  Re  =  1.25  x  10^  and  Re  =  1.2  x  10^,  are  plotted  in  Fig.  6  as  the 
solid  and  dashed  line  curves,  respectively.  Also  plotted,  for  reference,  is  a  straight 
line  with  a  slope  (logarithmic  derivative)  of  2/3 »  corresponding  to  the  high-Schmidt 
number  Batchelor  theoretical  spectrum.'* 


logiott 


Fig.  7  Frequency-scaled  concentration  spectra.  Solid  lines:  x/d  =  100  data  (1.25  < 
Re  X  10“*  <  7.?).  Dashed  lines:  x/d  =  305  data  (1.2  <  Re  x  10~*  <  6.5). 
Individual  spectra  are  offset  by  — 21ogio(Re/Reo),  Reo  =  1.2  x  10*,  for 
clarity. 


The  solid  lines  in  Fig.  7  plot  spectra  derived  from  measurements  at  x/d  =  100, 
for  Re  X  10“*  =  1.25,  1.76,  2.55,  3.6,  5.1,  and  7.2.  The  decrease  in  the  scaled 
frequency  resolved,  in  the  x/d  =  100  spectra,  as  the  speed  of  the  flow  increased 
with  increasing  Reynolds  number,  is  evident.  The  dashed  lines  plot  spectra  mea¬ 
sured  at  x/d  =  305,  for  Re  x  10“*  =  1.2,  2.4,  4.0,  and  6.5 .  As  can  be  seen,  the 
highest  frequency  resolved  in  the  x/d  =  305  spectra  is  a  much  weaker  function  of 
Reynolds  number,  signal-to-noise  ratio  limitations  being  more  serious  than  spatial 
resolution  at  this  station.  Individual  spectra  for  both  x/d  stations  are  plotted  offset 
by  — 21ogio(Re/Reo),  with  Rcq  =  1.2  x  10*,  to  aid  in  visualizing  the  evolution  of 
trends  with  Reynolds  number. 
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The  longer  characteristic  local  time  scale  at  xfd  =  305  did  not  allow  as  many 
large  scale  structures  to  be  included  in  the  record.  As  a  consequence,  the  statistical 
convergence  of  the  xfd  =  305  spectra  is  not  as  good  as  for  the  x/d  =  100  data. 
On  the  other  hand,  the  higlier  relative  spatial  I'esolution  at  xfd  —  305  allowed  the 
spectrum  to  be  estimated  to  a  higher  (scaled)  frequency.  That  trade-oif  aside,  the 
agreement  between  the  (scaled)  spectra  at  xjd—  100  and  x/d  =  305  holds  for  all 
the  cases  for  which  data  were  recorded  at  the  same,  or  nearly  the  same,  Reynolds 
number  at  the  two  stations. 


Fig.  8  Frequency-scaled  concentration  spectra  derived  horn  data  at  z/d  =  100  and 
1,25  <  Re  X  10"^  <  7.2  (no  offsets). 


Following  the  transition  out  of  the  large  scale  frequency  regime  (/tj  ~  1),  the 
spectra  appear  to  be  described  by  a  power-law  with  an  exponent  that  is  increasing 
from,  roughly,  1.2  to  1.5,  in  absolute  value,  with  increasing  Reynolds  number  (c/. 
Ref.  15,  Fig.  5.2).  This  progression  with  Reynolds  niunber  is  easier  to  discern  in 
Fig.  8,  which  plots  the  concentration  spectra  at  x/d  —  100,  for  Re  x  10“^  =  1.25 
(solid  line),  1.76  (dashed  line),  2.55  (dot-2-dash),  3.6  (2-dot-2-dash),  5.1  (dot- 
dash),  and  7.2  (2-dot-dash),  with  no  offsets.  The  extent  of  the  power-law  regime 
can  be  seen  to  increase  slightly  with  increasing  Reynolds  number.  A  similarly 
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increasing  spectral  slope  (in  absolute  value)  with  Reynolds  niunber  was  also  noted 
in  measurements  of  gas-phase  spectra,*^  for  Re  x  10~^  =  0.5,  1.6,  and  4.0. 

As  can  be  seen  by  sighting  along  the  spectra  in  Fig.  8,  the  power-law  region 
is  followed  by  a  different  regime  at  higher  frequencies  yet.  This  regime  does  not 
support  the  Batchelor  k~^  prediction^  that  should  apply  for  over  a  decade  and  a 
half  in  frequency  in  this  case  (recall  that  Sc  fas  1.9  x  10^  here).  This  can  be  seen  in 
Fig.  6,  which  includes  a  dotted  line  with  a  reference  slope  of  2/3  >  corresponding  to  a 
k~^  spectrum.  The  spectra  at  increasing  Reynolds  number  move  even  further  from 
this  slope,  as  can  be  seen  by  comparing  the  Re  =  1.2  x  10"*  spectnun  with  those  at 
higher  Reynolds  numbers  in  Fig.  8.  It  should  be  noted  that  this  conclusion  extends 
to  frequencies  below  frt  ~  2,  which  are  unaffected  by  resolution  and  compensation 
considerations. 


Fig.  9  Spectnun  slope  (logarithmic  derivative).  Solid  lines:  data  at  x/d  =  100  and 
Re  =  1.25  X  10*.  Dashed  line:  x/d  =  305,  Re  *=  1.2  x  10*. 


To  facilitate  the  study  of  this  higher  frequency  regime,  we  also  computed  the 
slope  of  the  spectra  (logarithmic  derivative).  A  plot  of  spectrum  slopes,  for  the 
lower  Reynolds  number  data,  appears  in  Fig.  9.  These  were  derived  from  the  two 
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x/d  =  100,  He  =  1.25  x  10^  spectra  in  Fig.  3  (solid  lines)  and  the  x/d  =  305, 
Re  =  1.2  X  10^  spectrum  (dashed  line).  Their  comparison  helps  assess  issues  of 
statistical  convergence  and  confidence  in  this  more  delicate  statistic,  as  well  as  the 
effects  of  spatial  resolution  and  the  applied  compensation. 

The  plots  in  Fig.  9  suggest  a  frequency-dependence  of  the  slope  of  the  spectra 
that  is  close  to  a  strzught  line,  in  these  coordinates,  for  frequencies  higher  thzm 
frg  ~  1.2,  at  this  Reynolds  number.  This  linear  beha\nor  extends  for  a  decade 
and  a  half  and  into  frequencies  beyond  which  the  data  are  limited  by  resolution  and 
signal-to-noise  considerations.  The  straight  line  appears  to  be  a  good  representation 
for  frequencies  below  frg  fa  2,  for  which  the  effects  of  compensation  were  negligible 
even  for  the  x/d  =  100  data  (c/.  Fig.  3).  For  frequencies  above  frg  ca  2,  the  same 
straight  line  also  describes  the  behavior  for  both  the  x/d  =  100  and  x/d  —  305 
data,  which  were  affected  by  resolution  (and  compensation)  to  a  different  extent 
(c/.  Figs.  3  and  5). 


logjot^Tj) 

Fig.  10  Slope  (logarithmic  derivative)  of  spectra  from  data  at  x/d  —  100  and  1.25  < 
He  X  10“*  <  7.2.  Spectra  offset  as  in  Fig.  7.  Line  types  as  in  Fig.  8. 
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A  straight  line  for  the  logarithmic  derivative  of  the  spectrum  corresponds  to  a 
spectrum  that  is  parabolic  in  log-log  coordinates,  or  log-normal  in  linear  coordinates, 
».c., 

Scifrs)  oc  exp|-i  [aln(/ri)  +  6]^|  .  (18) 

This  expression,  rather  than  a  k~^  power  law,  seems  to  be  the  appropriate  de¬ 
scription  of  our  jet  fluid  concentration  spectra  at  high  frequencies,  over  a  range  of 
frequencies  at  least  as  large  as  ,  *•£.,  a  decade  and  a  half,  in  this  case. 

Figure  10  plots  the  local  slope  (logarithmic  deri\'ative)  of  the  spectra  in  Fig.  8. 
The  offset  scheme  employed  in  Fig.  7,  and  line-types  employed  in  Fig.  8,  were  also 
used  here.  Straight  lines  can  be  seen  to  be  a  good  representation  for  the  spectrum 
slope  at  high  frequencies  with  Reynolds  number  dependent  values  of  the  param¬ 
eters  a  and  b  in  £q.  18.  The  end  of  the  pou'er-law  regime  zmd  the  beginning  of 
the  log-normal  range  can  be  seen  to  shift  to  higher  frequencies  with  increasing 
Reynolds  number.  Oxu:  data  admit  a  Kolmogorov  scaling  for  the  Reynolds  num¬ 
ber  dependence  of  this  transition  fiequency,  ».e.,  jRe*^^,  but  corresponding  to  the 
convection  (passage  frequency)  of  a  physical  scale  roughly  SO  times  larger  than  the 
estimated  Kolmogorov  scale  (computed  using  the  estimated  mean  centerline  energy 
dissipation^^  and  the  kinematic  viscosity).  This  transition  is  not  very  well  defined, 
however,  and  other  Reynolds  number  scaling  possibilities  cannot  be  ruled  out. 


5.  Conclusions 

This  work  has  investigated  temporal  scalar  (jet  fluid  concentration)  power  spec¬ 
tra  on  the  centerline  of  high  Schmidt  number  turbulent  jets,  in  the  Reynolds  number 
range  1.2  <  Re  x  10“^  <  7.2 .  Our  spectra  exhibit  a  power-law  regime  at  frequen¬ 
cies  above  the  local  large  scale  passage  frequency,  with  a  Reynolds  ntunber  depen¬ 
dent  exponent  increasing  (in  absolute  value)  from,  roughly,  —1.2  to  —1.5  over  the 
Reynolds  number  range  investigated.  This  corroborates  a  similar  finding  for  gas- 
phase  jet  fluid  concentration  spectra  measured  at  comparable  Reynolds  numbers.^^ 
At  higher  frequencies,  the  spectra  are  well  represented  by  a  log-normal  relation  with 
Reynolds  number  dependent  coefficients.  While  our  data  admit  a  Kolmogorov-like 


21 


scaling  for  the  beginning  of  the  log-normal  region  in  the  spectrum,  t.e., 
other  possibilities  cannot  be  ruled  out. 

We  appreciate  that  our  results  are  at  odds  with  the  classical  picture  of  high 
Schmidt  number  scalar  spectra.  We  do  not  find  a  -®/3  power-law  regime,  even 
though  our  measurements  were  conducted  at  Reynolds  numbers  where  such  behav¬ 
ior  has  previously  been  reported  for  high  Schmidt  number  jet  fluid  scalar  spectra. 
Finally,  despite  adequate  resolution  and  signal-to-noise  ratio,  our  data  do  not  sup¬ 
port  the  Batchelor  k~^  power-law  prediction.^  Specifically,  we  found  no  constant 
A;~^  slope  at  high  frequencies  and  a  spectral  slope  that  does  not  even  locally  attain 
a  value  of  —1. 

On  the  whole,  our  scalar  spectra  are  rather  similar  to  those  derived  by  Gargett 
from  ocean  measurements.^^  In  conjimction  with  her  data  and  analysis,  the  current 
results  raise  further  questions  about  the  universal  descriptions  of  scalar  spectra,  and 
their  applicability  to  some  of  the  canonical  flows,  such  as  turbulent  jets. 


Acknowledgments 

We  would  like  to  acknowledge  contributions  by  Dan  Lang  and  David  Dowling 
to  the  experiments.  This  work  was  supported  by  AFOSR  Grants  No.  83-0213  and 
88-0155,  and  GRI  Contract  No.  5087-260-1467. 


References 

^  Kolmogorov,  A.  N.,  “Local  Structure  of  IVirbulence  in  an  Incompressible  Viscous 
Fluid  at  Very  High  Reynolds  Nmnbers,”  Dokl.  Akad.  Nauk  SSSR  30,  299  (1941). 

^  Corrsin,  S.,  “On  the  spectrum  of  isotropic  temperature  fluctuations  in  isotropic 
turbulence,”  J.  Appl.  Pbys.  22,  469-473  (1951). 

^  Oboukhov,  A.  M.,  “Some  specific  featiues  of  atmospheric  turbulence,”  J.  Fluid 
Mecb.  13,  77-81  (1962). 


22 


*  Batchelor,  G.  K.,  “Small-scale  variation  of  convected  quantities  like  temperature 
in  turbulent  fluid.  Part  1.  General  discussion  and  the  case  of  small  conductivity,” 
J.  Fluid  Mech.  5,  113-133  (1959). 

®  Monin,  A.  S.,  and  Yaglom,  A.  M.,  Statistical  Fluid  Mechanics:  Mechanics  of 
Turbulence  II  (Ed.  J.  Lumley,  MIT  Press,  Cambridge,  MA,  1975). 

®  Gibson,  C.  H.,  “Kolmogorov  similarity  hypotheses  for  scalar  fields:  sampling 
intermittent  turbulent  mixing  in  the  ocean  and  galax}',”  Proc.  Roy.  Soc.  A  434, 
149-164  (1991). 

^  Chapman,  D.  R.,  “Computational  Aerodynamics  Development  and  Outlook,” 
AIAA  J.  17(12),  1293-1313  (1979). 

®  Batchelor,  G.  K.,  Howells,  I.  D.  ,  and  Townsend,  A.  A.,  “Small-scale  variation 
of  convected  quantities  like  temperature  in  turbulent  fluid.  Part  2.  The  case  of 
large  conductivity,”  J.  Fluid  Mech.  5,  134-139  (1959). 

*  Gibson,  C.  H.,  and  Schwarz,  W.  H.,  “The  tmiversal  equilibrium  spectra  of  tur¬ 
bulent  velocity  and  scalar  fields,”  J.  Fluid  Mech.  16,  365-384  (1963) 

Grant,  H.  L.,  Hughes,  B.  A.,  Vogel,  W.  M.,  and  Moilliet,  A.,  “The  spectrum  of 
temperature  fluctuations  in  tiu-bulent  flow,”  J.  Fluid  Mech.  34,  423-442  (1968). 

Gargett,  A.  E.,  “Evolution  of  scalar  spectra  with  the  decay  of  turbulence  in  a 
stratified  fiuid,”  J.  Fluid  Mech.  159,  379-407  (1985). 

Gibson,  C.  H.,  “Fossil  turbulence  and  intermittency  in  sampling  oceanic  mixing 
processes,”  J.  Geophys.  Res.  92,  C5,  5383-5404  (1987). 

Komori,  S.,  Kanzaki,  T.,  Murakami,  Y.,  and  Ueda,  H.,  “Simultaneous  measure¬ 
ments  of  instantaneous  concentrations  of  two  species  being  mixed  in  a  turbulent 
flow  by  using  a  combined  laser-induced  fluorescence  and  laser-scattering  tech¬ 
nique,”  Phys.  Fluids  A  1(2),  349-352  (1989). 

Miller,  P.  L.,  and  Dimotakis,  P.  E.,  “Reynolds  ntunber  dependence  of  scalar 
fluctuations  in  a  high  Schmidt  number  tiurbulent  jet,”  Phys.  Fluids  A  3(5),  1156- 
1163  (1991). 

Miller,  P.  L.,  Mixing  in  High  Schmidt  Number  Turbulent  Jets,  Ph.D.  thesis, 
California  Institute  of  Technology  (1991). 


23 


Dimotakis,  P.  E.,  and  Miller,  P.  L.,  ‘^Some  consequences  of  the  boundedness  of 
scalar  fluctuations,”  Phys.  Fluids  A  2(11),  1919-1920  (1990). 

Miller,  P.  L.,  and  Dimotakis,  P.  E.,  “Stochastic  geometric  properties  of  scalar 
interfaces  in  tinrbulent  jets,”  Phys.  Fluids  A  3(1),  168-177  (1991). 

Ware,  B.  R.,  Cyr,  D.,  Gorti,  S.,  and  Lanni,  F.,  “Electrophoretic  and  FVictional 
Properties  of  Particles  in  Complex  Media  Measured  by  Laser  Light  Scattering 
and  Fluorescence  Photobleaching  Recovery,”  Measurement  of  Suspended  Parti¬ 
cles  by  Quasi-Elastic  Light  Scattering  (Wiley,  NY),  255-289  (1983). 

Blake,  W.  K.,  Mechanics  of  Flow-Induced  Vibration.  H  (Academic  Press,  Or¬ 
lando,  FL,  1986). 

Yoimg,  F.  R.,  Cavitation  (McGraw-Hill,  London,  1989). 

Dowling,  D.  R.,  Lang,  D.  B.,  and  Dimotakis,  P.  E.,  “An  Improved  Laser-Rayleigh 
Scattering  Photodetection  System,”  Exp.  in  Fluids  7(7),  435-440  (1989). 

22  Chen,  C.  J.,  and  Rodi,  W.,  Vertical  Turbulent  Bue^’ant  Jets.  A  Review  of 
Experimental  Data  (Pergammon  Press,  Oxford,  1980). 

Dowling,  D.  R.,  and  Dimotakis,  P.  E.,  “Similarity  of  the  concentration  field  of 
gas-phase  turbulent  jets,”  J.  Fluid  Mech.  218,  109-141  (1990). 

2^  Dowling,  D.  R.,  Mixing  in  Gas  Phase  TlirbuJent  Jets,  Ph.D.  thesis,  California 
Institute  of  Technology  (1988). 

Friehe,  C.  A.,  Van  Atta,  C.  W.,  and  Gibson,  C.  H.,  “Jet  txn-bulence:  Dissipation 
rate  measurements  and  correlations,”  AGARD  Turbulent  Shear  Flows  CP-93, 
18.1-7  (1971). 

Becker,  H.  A.,  Hottel,  H.  C.,  and  Williams,  G.  C.,  “The  Nozzle-Fluid  Concen¬ 
tration  Field  of  the  Roxmd  Ibrbulent,  Free  Jet,”  J.  Fluid  Mech.  30(2),  285-303 
(1967). 


APPENDIX  B 


Dimotakis,  P.  E.  1993  “Some  issues  on  turbulent  mixing  and  turbulence,”  W. 
C.  Reynolds’  60‘**  birthday  TVirbulence  Symposium  (22-23  March  1993,  Monterey, 
CA),  GALCIT  Report  FM93-la. 


Some  issues  on  turbulent  mixing  and  turbulence 


Paul  E.  Dimotakis 

Graduait  Aeronautical  Laboratories 
California  Institute  of  Technology 
Pasadena,  California  911t5 


Abstract 

Recent  data  on  turbulent  mixing  suggest  that  the  mixing  transition,  previ¬ 
ously  documented  to  occur  in  shear  layers,  also  occurs  in  jets,  as  well  as  many 
other  flows,  and  can  be  regarded  as  a  universal  phenomenon.  The  resulting,  fully- 
developed  turbulent  flow  requires  a  minimum  Reynolds  number  of  Re  «  10*,  or  a 
Thylor  Reynolds  number  of  Rex  «  10*,  to  be  sustained.  Turbulent  mixing  in  this 
fully-developed  state  does  not  appear  to  be  universal,  however,  with  a  qualitatively 
different  behavior  between  shear  layers  and  jets. 
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1.  Introduction 

A  correct  description  of  turbulent  mixing  is  particularly  taxing  on  our  vmder- 
standing  of  turbulence;  such  a  description  relies  on  an  account  of  the  full  spectnun 
of  scales.  Specifically,  to  describe  the  entrainment  stage  that  is  responsible  for 
the  engulfinent  of  large  pockets  of  irrotational  fluid  species  into  the  turbulent  flow 
region,^  the  large  sceJe  flow  structures  need  to  be  correctly  described.  Secondly,  to 
describe  the  subsequent  kinematic  stirring  process  that  is  responsible  for  the  large 
interfacial  siuface  generation  between  the  mixing  species,  the  intermediate  range  of 
scales  must  be  correctly  accounted  for.  These  are  below  the  largest  in  the  flow  in 
size,  but  above  the  smallest  affected  by  viscosity  auid  molecular  diffusivity.  Finally, 
the  dynamics  at  the  smallest  scales  mtist  be  captured  to  describe  the  molecular 
mixing  process  itself.  These  three  phases  of  turbulent  mixing  were  clezirly  identifi  3d 
as  “more  or  less  distinct  stages”  in  the  1948  description  of  mixing  by  Eckart,^  who 
dubbed  them  as  the  initial,  intermediate,  and  final  stages,  respectively.  In  the  case 
of  mixing  of  high  Schmidt  munber  fluids,  i.e.,  fluids  cheiracterized  by  a  molecular 
diffusivity,  V,  that  is  much  smaller  than  the  kinematic  viscosity,  i/,  it  is  also  useful 
to  distinguish  between  the  vorticity-diffusion  stage,  whereby  velocity  gradients  are 
removed,  and  the  species-diffusion  stage,  which  removes  scalar  gradients.^ 

On  the  other  side,  successful  descriptions  emd  models  of  mixing  provide  us 
with  tests  of  aspects  of  turbulent  flow  that  are  difficult  to  probe  by  other  means 
—  experimentally,  nmnerically,  or  theoretically  —  at  the  high  Reynolds  munbers  of 
interest  here. 

An  exciting  discussion  in  the  context  of  mixing  by  chaotic  eidvection  has  been 
taking  place  during  the  last  decade.  See,  for  example.  Refs.  4-7,  as  well  as  several 
papers  from  the  1990  lUTAM  Symposium  on  Fluid  Mechanics  of  Stirring  and  Mixing 
(Published  in  Phys.  Fluids  A  5,  Part  2,  May  1991).  There  is  little  doubt  that  this 
progress  will  contribute  to  our  imderstanding  in  the  context  of  the  broader  issues 
of  tturbulent  mixing.  The  present  discussion,  however,  will  be  limited  to  flows  at 
Reynolds  numbers  that  are  high  enough  for  the  turbulence  to  be  regarded  as  fully- 
developed.  In  that  regime,  the  impact  of  the  recent  progress  in  the-  behavior  of 
deterministic,  chaotic  systems  has  yet  to  be  felt,  in  my  opinion. 
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As  a  practical  matter,  fully-developed  turbulent  flow  typically  requires  that  the 
local  flow  Reynolds  nmnber,  t.e., 

^  V(,x)Six) 

- Z — 

must  be  high  enough.  Generally  speaking,  txurbulence  cannot  be  sustained  if  the 
(local)  Reynolds  number  falls  below  some  minimum  value,  Remin-  the  expression 
in  Eq.  1,  the  characteristic  velocity,  U{x),  and  transverse  extent  of  the  flow,  6(x), 
are  to  be  tzdcen  as  local  values. 

The  main  part  of  the  discussion  will  be  drawn  from  turbulent  mixing  in  the  far 
fleld  of  nonaccelerating  (negligible  streamwise  pressure  gradient)  shear  layers  and 
jets.  While  these  two  flows  are  similar,  in  many  ways,  they  are  suflBciently  different 
in  others  to  be  useful  as  test  beds  of  ideas  and  prospects  for  universal  descrip¬ 
tions  of  ttirbulent  mixing  behavior.  The  phenomena  are  foimd  to  have  a  broader 
manifestation,  however,  with  conclusions  relevant  to  turbulent  flow  in  general. 

In  the  case  of  shear  layers,  the  characteristic  velocity  will  be  taken  as  the 
(constant)  freestream  velocity  difference,  t.c.,  f^sKa^)  =  fii(ac),  whereas  the 

characteristic  length  will  be  taken  as  the  local  shear  layer  width,  t.e.,  6si(z)  cx  x. 
Assuming  constant  fluid  properties,  t.e.,  v  —  const.,  this  yields  a  local  Reynolds 
number  for  shear  layers  that  increases  linearly  with  the  streamwise  coordinate,  t.e., 

Jie,i(x)  oc  X  .  (2a) 

In  the  case  of  jets,  the  charsM:teristic  velocity  will  be  taken  as  the  local  centerline 
velocity  of  the  jet,  t.e.,  Uy{x)  =  Uc(x)  oc  x~^,  while  the  local  length  scale  will  be 
taken  as  the  local  jet  diameter,  t.c.,  5(x)  =  ^j(x)  oc  x,  yielding  a  local  jet  Reynolds 
number  that  is  a  constant  of  the  flow,  t.e., 

Rej(x)  5^  fn(x)  .  (2b) 

This  difference  between  shear  layers  and  jets  in  the  dependence  of  the  local  Reynolds 
number  is  interesting  in  the  context  of  spatially  developing  flows  and  the  evolution 
of  the  distribution  of  scales  and  turbulence  spectra. 
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As  we  increase  the  flow  Reynolds  niimber,  from  small  values  to  values  approach¬ 
ing  some  minimum  Reynolds  number  for  fully-developed  txirbulence,  the  flow  is  able 
to  generate  ever-increasing  interfaciaJ  area  between  the  mixing  species,  increasing 
mixing  and,  in  the  case  of  chemically-reacting  flow,  chemical  product  formation  and 
heat  release.  Beyond  this  transition  region,  t.e.,  for  Re  ^  Remind  the  Reynolds  nrim- 
ber  dependence  of  the  amount  of  mixed  fluid  can  be  expected  to  be  weaker.  In  fact, 
a  tenet  of  fully-developed  turbulent  flow  theory  is  that,  at  high  enough  Reynolds 
numbers,  the  dependence  of  the  various  mean  flow  quantities  on  Reynolds  number 
should  become  negligible,  or  vanish. 

On  the  other  hand,  mixing  depends  on  the  behavior  of  gradients  in  the  flow  as 
well  as  concentration  of  diffusing  species  and  the  “principle  of  self-similarity  with 
respect  to  Reynolds  number  cannot  be  expected  to  be  apphcable  . . . ,  since  these 
gradients  are  determined  by  small-scale  fluctuations.”  (Ref.  8,  p.  xvi)  We  will  ex¬ 
amine  these  propositions  by  comparing  the  outcome  of  experimental  investigations 
on  turbulent  mixing  conducted  in  both  gas-  and  liquid-phase  shear  layers  and  jets. 


2.  Transition  Reynolds  numbers  in  shear  layers  and  jets 

A  qualitative  difference  in  the  appearance  of  the  scalar  field  is  observed  across 
the  transition  in  a  shear  layer  to  a  more  well-mixed  state,  as  the  Reynolds  number 
is  increased,  as  illustrated  in  Fig.  1.  The  transition  to  a  more  well-mixed  state  in 
a  gas-phase  turbulent  shear  layer  was  documented  by  Konrad,®  who  used  an  aspi¬ 
rating  probe^®  to  estimate  the  local  value  of  the  high-speed  stream  fluid  fraction, 
averaged  over  the  resolution  voliime  and  time-response  of  the  aspirating  probe.  Sub¬ 
sequent  estimates  of  mixing  and  chemical  product  volxune  fraction  in  liquid-phase 
shear  layers,**  using  a  pH  indicator,  as  well  as  estimates  derived  from  probability- 
density  fimctions  (pdf’s)  measmed  using  laser-induced  fluorescence  techniques,*^ 
documented  the  seune  behavior. 

The  results  from  the  two  liquid-phase  shear  layer  measurements  are  plotted 
in  Fig.  2,  which  depicts  the  estimated  chemical  product  thickness  as  a  fimction 
of  the  locaJ  Reynolds  number  at  the  measuring  station.  A  meurked  increase  in 
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Fig.  1  Laser-induced  fluorescence  streak  images  of  the  scalar  field  in  a  liquid-phase 
shear  layer,  for  Re  ~  1.75  x  10^  (left)  and  Re  ~  2.3  x  10^  (right).  Data  from 
Ref.  12,  Figs.  7  and  9,  respectively. 

the  estimated  chemical  product  can  be  seen  to  occur  at  Re  10^.  This  is  also 
associated  with  a  change  in  the  pdf  of  the  scalar  fluctuations.  In  the  pretransition 
region,  the  pdf  of  the  conserved  scalar  in  the  flow  is  dominated  by  the  near-delta- 
function  contributions  of  the  unmixed  (pure)  fluid  from  ezwih  of  the  freestreams.*^ 
In  the  posttrzmsition  regime,  the  composition  of  the  mixed  fluid  across  the  layer 
develops  a  preferred  value  that  is  well-correlated  with  the  one  inferred  from  the 
estimated  overall  entrainment  ratio  for  the  layer.^’^^  The  pdf  evolves  from  one  limit 
to  the  other  in  the  course  of  this  transition  (c/.  Ref.  12,  Sec.  5.4,  and  Ref.  13),  with 
a  relatively  long  memory  of  the  (typically,  much  larger)  initial  asymmetry  in  the 
relative  amounts  of  each  of  the  freestream  fluids. 

As  noted  in  the  discussions  of  these  experiments, finite  resolution  limita¬ 
tions  in  these  liquid-phase  experiments  overestimated  the  absolute  amount  of  chem¬ 
ical  product  by,  roughly,  a  factor  of  two,  as  confirmed  in  chemically-reacting  exper¬ 
iments  which  measured  the  chemical  product  volume  fraction  directly.*^  Neverthe¬ 
less,  the  documented  increase  in  the  iimount  of  mixing  at  the  transition  Reynolds 
munber  is  qualitatively  correct  auid  was  found  to  occur  at  the  same  Reynolds  number 
in  both  gas-  and  liquid-phase  shear  layers. 
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Fig.  2  Reynolds  ntunber  dependence  of  chemical  product  volume  fraction,  in  a 
liquid-phase  shear  layer,  in  the  vicinity  of  the  mixing  transition.*^  Note  that 
absolute  values  are  overestimated  in  both  experiments  (see  text). 

The  transition  to  a  more  well-mixed  state,  in  these  experiments,  was  correlated 
with  the  appearance  of  streamwise  vortices  and  the  ensuing  transition  to  three- 
dimensionality  of  flow  that  is  nominally  two-dimensional  in  the  initisd  region.**~*^ 
See  discussion  in  the  review  paper  by  Roshko.*^  Corroborating  evidence  was  also 
found  in  the  numerical  simulations  of  time-developing  shear  layers  by  Moser  & 
Rogers  that  followed  the  developing  flow  to  sufficiently  high  Reynolds  numbers  to 
document  the  beginning  of  this  behavior.*^ 

The  transition  to  a  more  well-mixed  state  in  turbulent  jets  is  less  conspicuous 
than  in  shear  layers.  TVirbulent  jets  being  three-dimensional,  even  at  low  Reynolds 
numbers,  such  a  transition  is  not  correlated  with  a  transition  to  three-dimensionality 
in  the  flow  fleld.  Nevertheless,  there  is,  again,  a  qualitative  difference  in  the  ap¬ 
pearance  of  the  scalar  fleld  for  values  of  the  Re3molds  number  that  are  lower  than 
Rcmin  ^  values  that  are  comparable  to  that,  or  higher.  This  is  illustrated 

in  the  laser-induced  fluorescence  images  in  Fig.  3,  of  the  jet-fluid  concentration  in 
the  plane  of  symmetry  of  liquid-phase  jets.**  Unmixed  reservoir  fluid  (black)  can  be 
seen  throughout  the  turbulent  region  and,  in  particular,  all  the  way  to  the  jet  axis 
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in  the  lower  Reynolds  number  (left)  image  at  Rej  ~  2.5  x  10^  (imaged  field  spans 
0  <  ar/d  <  35,  where  d  is  the  jet  nozzle  diameter).  This  is  not  the  case  in  the  higher 
Reynolds  number  (right)  image  at  Rej  ~  10^  (imaged  field  spans  0  <  x/d  <  200), 
in  which  jet  fiuid  of  varying  concentrations  can  be  seen  to  be  more  volume-filling 
within  the  turbulent  region. 


Fig.  3  Jet-fiuid  concentration  in  the  plane  of  symmetry  of  a  roimd  turbulent  jet. 
Left  image:  Rej  ~  2.5  x  10^  (0  <  x/d  <  35).  Right  image:  Rej  ~  10^ 
(0  <  x/d  <  200).  Data  from  Ref.  18,  Figs.  5  and  9. 


Seitzman  et  al}^  investigated  the  outer  entrainment  and  mixing  region,  using 
laser-induced  fiuorescence  images  of  OH  radicals  in  a  H2-air  turbulent  diffusion 
fiame.  A  qualitative  evolution  in  the  complexity  of  the  thin  btuning  regions  can 
be  seen,  as  the  Re3Tiolds  number  was  increased  from  2.3  x  10^  to  4.95  x  10^  (c/. 
their  Fig.  3).  In  these  experiments,  this  evolution  is  also  influenced  by  decreasing 
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Fig.  4  Normalized  variance  of  jet-fluid  concentration  on  the  axis  of  a  turbulent 
jet,  as  a  function  of  the  jet  Reynolds  number  (Ref.  22,  Fig.  7.2).  Circles: 
Liquid-phase  jets.^^  TViangles:  Gas-phase  jets.*° 

buoyancy  and  the  decreasing  relative  importance  of  baroclinic  vorticity  generation, 
as  the  Reynolds  number  was  increased,  and  is,  therefore,  not  entirely  attributable 
to  Reynolds  number  effects. 

Similar  behavior  is  reflected  in  the  measurements  of  the  rms  of  the  scalar  (jet 
fluid  concentration)  fluctuations  on  the  axis,  in  the  far  field  of  gas-  and  hquid- 
phase  jets,  as  a  function  of  jet  Reynolds  number.^®"^*  The  data,  in  the  form  of  the 
normalized  scalar  fluctuation  variBmce,  are  plotted  in  Fig.  4  (Ref.  22,  Fig.  7.2).  The 
liquid-phase  data  exhibit  a  decrease  in  the  fluctuation  level  with  Reynolds  number, 
with  a  rather  less  sensitive  dependence  for  Reynolds  numbers  higher  than  Re 
2  X 10^,  or  so.  Noting  that  lower  fluctuation  levels  correspond  to  more  homogeneous 
mixing,  t.e.,  a  pdf  of  concentration  values  that  are  more  tightly  clustered  around 
the  local  mean,  we  see  that,  at  least  for  the  case  of  a  liquid-phase  jet,  the  flow 
transitions  to  a  more  well-mixed  state  as  the  Reynolds  number  is  increased,  as  in 
the  shear  layer,  even  though  in  a  more  gradual  manner  (c/.  Fig.  2).  A  much  weaker 
Reynolds  niunber  dependence  of  the  normalized  scalar  variance  can  be  seen  for  the 
gas-phase-jet  data. 
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The  Reynolds  number  dependence  of  turbulent  mixing  and  chemical  product 
formation  in  turbulent  jets  was  recently  investigated  in  gas-phase  jets.^^'^^  In  this 
context,  the  turbulent  diffusion  flame  length,  Z/f,  is  important  and  marks  the  dis¬ 
tance  from  the  nozzle  required  to  mix  and  bum  the  reactant  carried  by  the  jet 
fluid.  If  the  stoichiometry  of  the  jet/reservoir  reactants  and  jet  entrainment  are 
held  constant,  and  for  fast  chemical  kinetics  (high  Dmnkohler  number  limit),  the 
flame  length  dependence  on  the  various  flow  parameters  provides  us  with  a  measure 
of  the  dependence  of  mixing  on  those  parzuneters.  Decreasing  flame  lengths,  for 
example,  imply  faster  (better)  mixing. 

The  dependence  of  the  flame  length  on  the  stoichiometry  of  the  jet-/reservoir- 
fluid  chemical  system  must  first  be  factored  in  the  analysis.  In  particular,  for  a 
momentum-dominated,  turbulent  jet  diffusion  flame,  the  flame  length  is  linearly 
dependent  on  the  (mass)  stoichiometric  mixture  ratio  (e.g.,  Refs.  25-27),  i.c., 

^  ~  -h  S  ,  (3) 

where  is  the  mass  of  reservoir  fluid  required  to  completely  consume  a  \mit  mass  of 
jet  fluid.  The  measturements  must  then  be  regarded  as  investigations  of  the  behavior 
of  the  stoichiometric  coefficient.  A,  and  normalized  virtual  origin  (intercept),  B,  and 
their  dependence,  in  turn,  on  the  flow  parameters. 

In  these  experiments,  long  platinum  wires  were  stretched  across  the  tiurbulent 
diffusion  flame  and  spaced  in  equal  logarithmic  increments  along  the  jet  axis.  These 
permitted  the  line-integral  of  the  temperatinre  rise,  AT(x,  y),  due  to  heat  released 
in  the  chemical  reaction,  to  be  measured  along  the  y-coordinate  (transverse  to  the 
jet  axis),  as  a  function  of  the  downstream  coordinate.  See  Fig.  5. 

The  experiments  utilized  the  F2  -l-NO  chemical  reaction,  with  F2  diluted  in  N2 
forming  the  jet  fluid,  and  NO  diluted  in  N2  forming  the  quiescent  reservoir  fluid. 
With  this  chemical  system,  an  adiabatic  flame  temperatmre  rise,  ATf,  as  low  as  7K 
was  realized  (with  the  reaction  still  in  the  fast-kinetic  regime).  Such  low  values 
were  dictated  by  the  results  of  a  separate  investigation  that  assessed  the  effects  of 
buoyancy  and  ascertained  that  the  measurements  were  realized  in  the  momentum- 
dominated  regime  for  this  heat-releasing  flow.  The  Reynolds  number  was  varied  by 
varying  the  pressure  in  the  combustion  vessel. 
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stoichiometric  mixture  ratio,  The  data  plot  the  product  thickness  Sp,  normal¬ 
ized  with  the  length,  I>w  i  of  the  platinum  resistance  wire  used  to  measure  the  line 
integral,  versus  the  logarithm  o{  x/d*,  where  d*  is  the  jet  source  diameter.^^ 


1.4  1.6  1.8  2.0  2.2  2.4  2.6 

lOOioix/a*) 

Fig.  6  Product  thickness  vs.  log|o(a:/d*),  for  several  stoichiometric  mixture  ratios 
(Ref.  24,  Fig.  5). 

To  analyze  these  data,  we  note  that  for  regions  of  the  flow  well  upstream  of 
the  flame  length,  t.e.,  for  x  ^  Li,  the  entrained  reactant  is  consumed  on,  or  just 
outside,  the  boimdary  of  the  turbulent  region.  There,  the  tmrbulent  fluid  is  jet- 
fluid-reactant  rich  and  it  need  comprise  only  a  small  fraction  of  the  mixed  fluid  to 
consume  the  entrained  reservoir-fluid-reactant.  The  diffusion/reaction  process  then 
takes  place  in  a  relatively  thin  peripheral  reaction  zone  at  y  =  ±ilr(a;)i*  whose 
ensemble-averaged  radius,  Rt{x),  is  proportional  to  x.  As  a  consequence,  the  line 
integral  of  the  time-averaged  temperature  rise  across  the  tinbulent  region  increases 
as  the  chemical  reaction  releases  heat  in  the  thin  reaction  zones  at  the  edges  of  the 
turbulent  region. 

*  This  picture  is  corroborated  by  the  OH-imsges  obtained  by  a  number  of  investigators  in  H2-air 
jet  flames  (c/.  Refs.  19  and  31,  for  example). 
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It  was  conjectured  that  the  radial  integral  of  the  temperature  rise,  at  a  given 
station  x,  increases  in  proportion  to  the  entrainment  velocity  at  that  station, 

tie[iJr(x)],  i.C., 


^  /  AT(l,r)dr  <x  , 

or,  for  a  momentum-driven,  turbulent  jet. 


/“ 

OO 


Ar(x,y) 


dy  a 


dx 


dx 


QC 


-oo  ATf  ■■  27ri2,(x) 

Integrating  this  relation  and  scaling  with  the  flame  length  £f,  we  have 

^p(x  ^  L{) 


^  (i) 


+  b  , 


(5) 

(6) 

(7) 


This  dependence  of  the  hne  integral  on  x  suggested  the  logarithmic  wire  spacing 
used  in  the  experiment  and  was  used  in  the  analytical  form  of  the  flt  for  the  line- 
integrated,  time-averaged,  temperature-rise  data  (c/.  Fig.  6). 


Beyond  the  end  of  the  flame  region,  t.e.,  for  x  >  Lt,  no  further  heat  is  released 
and,  in  the  absence  of  buoyancy  effects,  the  temperature  excess  becomes  a  passively- 
convected  scalar  with  a  self-similar  proflle.  In  that  case, 

Ar(,,»)  =.  AT(*.O)/0  CX  i/0 

and  the  product  thickness  hne  integral  becomes  independent  of  the  downstream 
coordinate,  x,  i.e., 

Sp(x  >  Lt)  oc  f 

•/— OO 


i  /  (j)  d»  ^  b(x)  .  (8) 


As  can  be  seen  in  Fig.  6,  the  experimental  results  conflrmed  the  conjecture  for 
X  Li.  They  were  also  consistent  with  the  anticipated  conserved-scalar  behavior 
of  the  temperature  rise  for  x  >  Lf,  i.e.,  a  product  thickness  that  asymptotes  to  a 
constant  value.  Such  data  allow  us  to  estimate  the  flame  length,  Lt.  In  particular, 
one  ca.n  accept  an  operational  definition  of  Lt  as  the  location  where  the  product 
thickness  hne  integral  (Ek).  4)  has  attained  99%  of  its  asymptotic  value,,  as  one  does 
on  the  basis  of  a  boimdary  layer  velocity  profile,  for  example. 
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Figure  7  (Ref.  23,  Fig.  4.8)  plots  the  stoichiometric  coefl5cient,  A,  in  the  ex¬ 
pression  for  the  flame  length  (c/.  Eq.3),  i.c.,  the  slope  of  the  flame  length  vs.  the 
stoichiometric  mixtme  ratio  This  can  be  regarded  as  the  additional  length,  in 
units  of  the  jet  somce  diameter  d*,  required  to  entrain,  mix,  and  react  with  a  unit 
increase  in  the  stoichiometric  ratio  of  the  jet-/reservoir-fluid  chemical  system,  i.e.. 


In  the  fast  kinetic  regime,  as  was  the  case  in  these  expenments,  this  quantity  is  a 
useful  measme  of  mixing.  It  separates  the  self-similar,  far-field  behavior  from  that 
of  the  virtual  origin  in  the  overall  mixing  process. 


•f  - - 

3.5  4.0  4.5  5.0  5.5 

log^QRe 


Fig.  7  Flame  length  stoichiometric  coeflicient  i4  (Eq.  3).  Squares:  Geis- phase  chemi¬ 
cally-reacting  jet  data.^^  Diamond:  Laser-induced  fluorescence,  liquid-phase 
data.®^  Lambda:  pH-indicator,  liquid-phase  data.^®’^^  Triangles:  Flame 
length  data  inferred  from  gas-phase,  nonreacting  data  (see  text).^^ 

The  data  in  Fig.  7  indicate  that  mixing  in  the  far  field  of  a  tmbulent  jet  improves 
relatively  rapidly  with  increasing  Reynolds  number.  Specifically,  A  decreases  until  a 
Reynolds  number  of,  roughly,  2  x  10^,  with  a  much  weaker  dependence  on  Reynolds 
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number,  if  any,  beyond  that.  These  data  are  in  accord  with  the  nonreacting,  liquid- 
phase  data  in  Fig.  4,  which  also  indicate  improved  mixing  up  to  Reynolds  numbers 
of  2  X  10^,  or  so,  with  a  weaker  dependence  beyond  that.  The  latter  data,  however, 
do  not  permit  the  separation  of  the  far-field  and  virtual-origin  contributions  to  the 
overall  mixing  process,  as  do  the  chemically  reacting  jet  data.  We  should  also  note 
that  the  near-  and  intermediate-field  behavior,  which  contributes  to  the  virtual 
origin  of  the  mixing  process  and  the  resulting  fiame  length,  does  not  exhibit  the 
same  Reynolds  number  dependence.^* 

A  potential  difficulty  should  be  recognized  between  the  inferred  behavior  based 
on  the  nonreacting,  gas-phase  jet  data  (triangles),*^  and  the  chemically-reacting, 
gas-phase  data  (squares).**  Two  observations  are  relevant  here.  The  values  es¬ 
timated  from  the  nonreacting  gas-phase  data  (triangles)  were  derived  assuming 
certain  similarity  properties  of  the  concentration  pdf  zind  the  virtual  origin  of  the 
whole  process  (see  discussion  in  Ref.  34,  Sec.  5.4,  and  Ref.  21,  Appendix  B). 

Partly  as  a  consequence,  as  was  also  noted  in  the  comparison  between  the  data 
in  Figs.  4  and  7,  it  is  not  possible  to  separate  the  contribution  to  the  flame  length 
of  the  (rather  large)  virtual  origin  of  the  mixing  process,  and  its  dependence  on 
Reynolds  niunber,**  from  the  Reynolds  number  dependence  of  the  far-field  mixing 
process,  ».c.,  the  dependence  of  the  flame  length  stoichiometric  coefficient,  A. 

While  this  may  be  a  minor  point,  we  may  wish  to  note  that  transition  Reynolds 
numbers  for  jets  seem  to  be  twice  as  large  as  for  shear  layers.  On  the  one  hand, 
the  two  flows  are  sufficiently  different  to  admit  differences  in  their  behavior  of  a 
factor  of  2,  or  so,  in  Reynolds  number.  On  the  other,  however,  if  the  characteristic 
large  scale  6{x)  chosen  for  the  local  Reynolds  number  definition  of  a  jet  is  the  local 
radius ,  as  would  be  appropriate  if  the  length  scale  in  the  general  case  is  defined  as 
the  treinsverse  spatial  extent  across  which  the  shear  is  sustained,  then  the  transition 
Reynolds  number  for  jets  becomes  very  close  to  that  for  shear  layers. 
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3.  Mixing  in  fully-developed  turbulence 


In  fully-developed  turbulent  shear  layers,  beyond  the  mixing  transition,  mixing 
also  exhibits  a  weedcer  dependence  on  Reynolds  number.  Figures  8a  and  8b  plot 
experimentally  obtained  data  of  the  normalized  chemical  product  thickness,  t.e.. 


^  J_  p  Ar(x,y;^) 


(10) 


(c/.  Eq.4),  for  both  low  (Fig.  8a)  and  high  (Fig.  8b)  values  of  the  stoichiometric 
mixture  ratio  <f>.  The  stoichiometric  mixttire  ratio,  in  this  case,  denotes  the  parts 
(moles)  of  high-speed  fluid  required  to  consume  a  part  (mole)  of  low-speed  fluid  in 
the  shear  layer  mixing  zone  (e.g  .,  Ref.  30).  The  product  thickness  data  in  Figs.  8a 
and  8b  were  normalized  by  the  1%  temperature  rise  product  thickness,  67,  which 
has  been  found  to  be  close  to  the  visual  thickness,  6vis>  of  the  shezu-  layer.^’t^’^^*’^® 


4.0  4.5  5.0  5.5  6.0  6.5  7.0 

log  (Rej.^.) 

Fig.  8a  Low-^  normalized  product  thickness  w.  Reynolds  number  for  a  turbulent 
shear  layer.  Matched  firee-stream  density,  incompressible,  gas-phase  shear 
layers:  <i>  =  1/S  (square),  ^  =  1/4  (circle).^®  VerticsJ  lines:  <f>  =  1/8.’® 
Triangles:  </>  =  1/8.’^  Supersonic  shear  layers:  4>  =  1/4,  Pa/pi  =  0.71, 
Me  =  0.51  (arrow);  <f>  =  1/3,  P2/P1  =  5.95,  Me  =  0.96  (^amond).’* 
Liquid-phase  shear  layer:  4>  —  1/10  (cross). 
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log  (Reg^) 

Fig.  8b  High-^  shear  layer  normalized  product  thickness  vs.  Reynolds  number. 
Matched  free-stream  density,  incompressible,  gas-phase  shear  layers:  ^  =  8 
(square),  <f>  =  4  (circle).^®  THangles:  <l>  =  8.®^  Supersonic  shear  layers: 
^  =  4,  P2/P1  =  0.71,  Me  =  0.51  (arrow);  <^  =  3,  P2/P1  =  5.95,  Me  =  0.96 
(diamond).®*  Liquid-phase  shear  layer:  ^  =  10  (crosses).*® 

These  data  were  all  derived  from  chemically-reacting  shear  layer  experiments, 
conducted  in  the  fast  chemical-kinetic  regime,  t.e.,  high  Damkohler  munber  limit. 
In  this  limit,  all  molecularly  mixed  fluid  produces  chemiceJ  product,  as  dictated  by 
the  stoichiometry  of  the  mixed  fluid  composition  (c/.  discussion  in  Ref.  30,  Sec.  IV). 


The  chemically  reacting  experiments  were  conducted  for  low  and  high  values 
of  the  (molar)  stoichiometric  mixtxure  ratio,  that  can  be  expected  to  yield  near- 
stationary  values  of  the  product  thickness  ^p(^),  with  respect  to  <f>.  The  chemical 
product  thickness  has  a  rather  sensitive'  dependence  on  as  the  stoichiometric 
mixture  fraction,  =  <f>/{<f>+l)y  approaches  zero  or  miity,  as  a  result  of  the  regions 
near  the  delta  functions  correi^ponding  to  the  pure  fluid  in  the  composition  pdf.®®’*° 
Data  points  indicated  by  crosses  were  derived  from  experiments  in  liquid-phase 
flows.*®  All  other  data  were  based  on  results  from  gas-phase  flows.  The  plotted 
points  based  on  the  experiments  by  FVieler®®  and  Hall  ei  of.®*  were  calculated  by 
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Chris  Bond  from  the  original  temperature-rise  data.**  Data  at  the  highest  values 
of  the  Reynolds  number  were  derived  from  chemicsdly-reacting,  supersonic  sheeir 
layers  (Mj  >  1,  M2  <  1)^®  and,  as  a  consequence,  c£in  be  eissumed  to  be  susceptible 
to  the  combined  effects  of  compressibility  as  well  as  Reynolds  number.  We  note 
here  that  even  though  the  velocity  ratio,  r  =  U2/U1,  for  the  supersonic  shear  layer 
experiments  was  very  low,  and  should  perhaps  not  be  included  in  the  same  plot 
as  the  subsonic  shear  layer  experiments  for  which  r  ~  0.4,  the  molar  entrainment 
ratio,  Ea,  t.e.,  the  ratio  of  the  number  of  high-speed  freestream  fluid  moles  per 
mole  of  low-speed  freestream  fluid  entrained,  was  estimated  to  be  close  to  its  value 
for  the  incompressible  flow  experiments.  Specifically,  for  the  incompressible  flow 
experiments,®’^^’®®  En  ~  1.3.  The  estimated  values  for  the  compressible  shear 
layers  were  E^  1.07,  for  the  Me  0.51  shear  layer  (r  ~  0.24),  and  Ej^  ~  1.2,  for 
the  Me  ^  0.96  shear  layer  (r  ~  0.10).®*’^®’^® 

Several  observations  can  be  made  on  the  basis  of  these  data: 

a.  In  the  limit  of  fast  kinetics,  the  chemical  product  formed  decreases  as 
the  Reynolds  number  increases.  This  implies  less  efficient  mixing  in  fully- 
developed  shear  layers  as  the  Reynolds  number  increases.  Presently  avail¬ 
able  data  suggest  that  this  is  opposite  the  behavior  exhibited  by  turbulent 
jets  (recall  data  in  Figs.  4  and  7,  and  related  discussion). 

b.  The  data  admit  a  Reynolds  number  dependence  of  the  chemical  product 
thickness  for  the  low-^  case  (Fig.  8a)  that  is  stronger  than  the  high-<^  case 
(Fig.  8a). If  this  is  proven  to  be  the  case,  it  would  imply  a  complicated 
dependence  of  the  pdf  of  compositions  on  Reynolds  number. 

c.  Liquid-phase  sheau'  layer  mixing  exhibits  a  weziker  Reynolds  number  de¬ 
pendence  than  gas- phase  shear  layer  mixing  (c/.  Fig.  8b).  Unless  the  data 
foi  turbulent  jets  should  all  be  regarded  as  not  in  ftiUy-developed  flow  (c/. 
Fig.  4),  this  is  also  opposite  the  behavior  foimd  in  turbulent  jets. 


Private  communication. 

^  Recall  that,  for  these  flows,  the  high-speed  stream  fluid  is  preferentially  entrained.® 
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Based  on  the  presently  available  scant  experimental  data  and  vinless  the  region 
2  X  10*  <  Re  <  10®  conceals  a  second  transition  to  a  Reynolds-niimber-independent 
turbulent  mixing  state, 

d.  the  trend  in  Figs.  8a  and  8b  suggests  that,  at  least  at  low  compressibility, 
e.g.,  up  to  the  lower  of  the  two  convective  Mach  numbers  investigated, 
Reynolds  number  effects  dominate  compressibility  effects. 

The  latter  observation  must  be  regarded  as  the  most  tentative  and  must  await  the 
results  of  further  experimental  investigations  to  be  confirmed. 

It  is  interesting  that  in  the  case  of  mixing  in  turbulent  jets,  gas-phase  data 
for  the  flame  length  indicate  only  a  weak  Reynolds  number  dependence  beyond  the 
mixing  transition,  if  any  (Fig.  7).  Measurements  in  liquid-phase  jets,  however,  at 
Reynolds  numbers  as  high  as  7.2  x  10^,  yield  scalar  spectra  that  have  not  converged 
to  a  Reynolds-number-independent  state,^*  in  accord  with  the  data  in  Fig.  4. 

Finally,  it  is  noted  that  some  of  these  observations  are  at  variance  with  the  in¬ 
ferences  drawn  in  a  review  by  Broadwell  &:  Mungal,^^  of  earlier  data.  The  interested 
reader  is  directed  to  that  discussion  for  further  details. 


4.  TVansition  Reynolds  numbers  in  other  flows 

The  observations  of  mixing  transitions  in  shear  layers  and  jets  suggest  that  a 
minimum  Reynolds  number  may  be  required  for  turbtilence  to  develop  into  a  more 
well-mixed  state  in  these  flows.  Specifically,  we  must  have  Re  >  Reoin)  with  Remin 
in  the  neighborhood  of  0.5  x  10^  to  2  x  10®,  for  fully-developed  turbulent  flow.  It  is 
interesting  that  this  value  does  not  appear  to  be  peculiar  to  the  far-field  behavior 
in  turbulent  jets  and  free-shear  layers.  Other  flows  also  exhibit  similsu*  transitions 
at  comparable  vjdues  of  the  Reynolds  number,  as  we’ll  discuss  below. 

Pipe  flow,  for  example,  transitions  out  of  its  slug/puff  regime  to  a  less  inter¬ 
mittent,  fuUy-turbulent  state  over  a  range  of  Reynolds  numbers  that  depend  on  the 
entrance  conditions.  This  sensitivity  to  initial  conditions  diminishes,  however,  at  a 
Reynolds  number  in  the  vicinity  of  10®  .®* 
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The  Coles’  turbiilent  boundary  layer  wake  parameter,  II,  that  scales  the  outer 
flow  region  of  a  turbulent  boundsry  layer, is  foimd  to  increase  with  Reynolds 
number,  for  a  zero-pressure-gradient  boundary  layer,  attaining  an  asymptotic  value 
of  n  =  0.620  at  a  Reynolds  number,  based  on  displeurement  thickness,^^  of  Re  = 
UqoS*  fv  ~  0.8  X  10^.  See  Refs.  48  zind  49,  for  a  discussion,  and  Ref.  49,  Table  4 
and  Fig.  6,  for  a  compilation  of  low-speed,  turbulent  boimdary-layer  flow  data. 

In  experiments  by  Liepmann  &  Gharib,  in  the  near  field  of  turbulent  jets,  the 
number  of  azimuthal  nodes  in  vortex  structures  becomes  difficult  to  identify  beyond 
a  certain  Reynolds  number,  where  the  flow  transitions  to  a  much  more  chaotic 
state.®®  The  authors  correlate  this  transition  with  a  laminar-turbulent  transition  in 
the  jet  nozzle  boimdary  layers.  It  is  also  interesting,  however,  that  it  occurs  at  a 
Reynolds  number  very  close  to  10^. 

In  recent  experiments  on  lifted-flame  behavior.  Hammer  notes  a  change  in 
the  sc2ded  lift-off  height  of  turbulent  jet  flames  at  a  jet  Reynolds  number,  in  the 
neighborhood  of  Re  «  1.8  x  10^,  beyond  which  the  Reynolds  nvunber  dependence 
is  weaker.  See  data  in  Ref.  51,  Fig.  3.8,  and  discussion  following. 

In  his  review  of  bluff-body  flows,  Roshko  dociunents  several  regimes,  as  indi¬ 
cated  by  the  behavior  of  the  base  pressure  of  a  circular  cylinder,  as  a  function  of 
Reynolds  number.®^  In  particular,  the  (negative)  base  pressrire  is  foimd  to  increase 
in  the  range  of  Reynolds  numbers  of  0.3  x  10^  <  Re  =  Uoo  dzyxfv  <  2  x  10^  (Ref.  52, 
Fig.  1).  Roshko  attributes  this  behavior  to  a  transition  in  the  separating  shear 
layers. 


Measurements  of  the  scaled  turbulent  kinetic  energy  dissipation  rate. 


a 


(11) 


in  flow  behind  square  grids,  where  e  is  the  kinetic  energy  dissipation  per  imit  mass, 
(.  is  the  longitudinal  length  scale,  and  u'  is  the  rms  strejunwise  velocity  fluctuation 
level,  suggest  that  it  decreases  relatively  rapidly  with  increasing  Taylor  Reynolds 
number. 
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where  At  is  the  Taylor  microscale,  until  a  value  of  Rej  fa  70,  and  then  becomes 
much  less  sensitive  to  Reynolds  munber,  attaining  a  value  of  a  »  1  at  higher 
Reynolds  niunbers.  This  value  may  not  be  universal,  however,  with  measured  values 
in  the  range  1  <  a  <  2.7  behind  non-square  grids.®* 

A  similsu’  conclusion  was  arrived  at  by  Jimenez  ei  in  their  nmnerical 
simulations  of  tmbulence  in  a  spatially-periodic  cube,  in  the  range  36  <  Rer  <  170. 
They  report  a  value  of  a  ~  0.65  attained  for  Rer  >  95.  Since 

Rer  «  Re^^^  ,  (12b) 

it,  again,  appears  that  Re  >  Remin  fa  10^  is  a  necesseiry  condition  for  fully-developed 
tmbulent  flow.*® 

In  thermal  convection,  a  transition  from  “soft  turbulence”  to  “hard  turbulence” 
has  been  noted  for  Rayleigh  numbers  given  by  Ra  «  10®,  that  is  marked  by  a 
qualitative  change  in  the  pdf  of  the  measured  temperature  fluctuations.®®  Since 
Re  «  Ra^/*  for  this  flow,®^  we  again  recover  a  minimum  Reynolds  number  boimdary 
of  the  fully-developed  turbulent  state  at  Re  «  10®. 

Careful  experiments  were  recently  performed  that  measmred  the  torque  in 
Couette-Taylor  flow,  in  the  range  of  Reynolds  numbers  of  800  <  Re  <  1.23  x 
10®  These  experiments  revealed  a  “well-defined,  uon-hysteretic  transition”  in 
a  narrow  range  of  Reynolds  numbers,  10®  <  Rctr  <  1.3  x  10®.  The  flow  was  foimd  to 
be  qualitatively  different,  below  and  above  this  transition,  as  illustrated  in  the  flow- 
visualization  data  reproduced  in  Fig.  9,  with  pre-  and  post-transition  differences 
reminisce-it  of  the  corresponding  ones  in  jets  (c/.  Fig.  3).  See  also  additional  flow- 
visualization  data  in  Ref.  58,  Figs.  la,b.  Beyond  this  transition,  the  dependence 
of  the  torque  on  Reynolds  number  becomes  progressively  weaker.  Significantly, 
however,  the  torque  does  not  attain  viscosity-independent  behavior  to  the  highest 
Reynolds  numbers  investigated. 
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Fig.  9  Couette- Taylor  flow-visualization  data  at  (a)  Re  =  0.6  x  10“*,  and  (b)  Re  = 
2.4  X 10^.  From  Ref.  59,  Figs.  5a, b,  reproduced  by  kind  permission  of  Prof.  H. 
Swinney. 

5.  A  criterion  for  fully-developed  turbulence? 

The  preceding  observations  suggest  there  may  exist  a  property  of  turbulence 
that  induces  it  to  transition  to  a  more  well-mixed  state,  is  associated  with  Reynolds 
numbers  in  excess  of  Rcmin  ^  10^ »  appears  to  be  rather  independent  of  the 
details  of  the  flow  geometry.  The  following  is  a  proposed  ansatz  to  account  for  this 
behavior. 

That  this  transition  appears  to  be  independent  of  the  flow  geometrj'  indicates 
that  the  explemation  should  not  be  sought  in  the  large-scale  dynamics,  or  the  de¬ 
velopment  of  distinct  features  zmd  organized  patterns  in  these  flows.  These  are, 
typically,  flow-geometry  dependent.  One  is  rather  led  to  consider  the  physical  sig¬ 
nificance  of  the  various  scales  of  turbulence  and  their  Reynolds  number  scaling. 
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The  1941  Kolmogorov  proposals  hypothesize  that  the  dynamics  in  the  (inertial) 
range  of  scales  A  that  are  unaffected  by  the  outer  scale  6,  but  are  large  compared 
to  the  inner,  dissipation  (Kolmogorov)  scale, 

Ak  -  (-)  ,  (13) 

t.e.,  for 

Ak  <  A  <  ^  ,  (14) 

can  be  treated  in  a  universal,  self-similar  fashion.  In  this  range  of  scales,  for  example, 
the  energy  spectrum  is  predicted  (and  foimd)  to  exhibit  a  power-law  behavior  with 
a  -  S/a  exponent.®® 

To  refine  the  bounds  in  Eq.  14,  we  appreciate  that  independence  from  the  dy¬ 
namics  of  the  outer  scale,  6,  requires  that  the  scale  A  be  smaller  than  a  scale  that 
can  be  generated  directly  from  the  outer  scale  6.  Such  a  scale  would  an  outer  lam¬ 
inar  layer  thickness,  Al,  that  can  be  generated  by  a  single  6-size  sweep  across  the 
whole  transverse  extent  of  the  turbulent  re^on,  for  example.  The  size  of  this  scale 
can  be  estimated  in  terms  of  the  99%  thickness  of  a  Blasius  boimdary  layer,  for 
example,  that  is  growing  over  a  spatial  extent  6,  *.e., 

^  .  (15) 

0 

It  is  a  scale  connected  by  viscosity  to  the  outer  scale,  6,  of  the  flow.  By  virtue  of  its 
dependence  on  Reynolds  number  (c/.  Eqs.  12  and  15),  this  scale,  as  noted  by  H.  W. 
Liepmann  in  private  conversation  many  years  ago,  is  closely  related  to  the  Taylor 
microscale.  At- 

At  the  other  end  of  the  spectnun,  the  requirement  that  the  motions  must  be 
inviscid  with  respect  to  the  inner  dissipation  scales  dictates  that  the  local  scale 
A  must  be  large  with  respect  to  an  inner  viscous  scale,  A„  (c/.  Fig.  10)  that  can 
be  taken  as  proportional  to  the  (defined)  Kolmogorov  dissipation  scale,  Ak-  This 
allows  us  to  refine  the  inequality  that  boimds  the  inertial  range  of  scales  (Eq.  14) 
to  the  one  below,  *.c., 


(16) 


Fig.  10  Schematic  of  the  outer  scale,  6;  the  Taylor  scale,  At;  and  the  viscous  scale, 
A„,  in  a  sheared  turbulent  region. 

as  a  necessary  condition  for  fully-developed  flow. 

To  translate  this  inequality  to  a  relation  for  the  Reynolds  number,  we  need  the 
Reynolds  number  dependence  of  the  ratio  of  the  various  scales  to  the  outer  scale 
6.  We  can  rely  on  Eq.  15  for  the  estimate  of  the  outer  laminar-layer  thickness,  Al, 
suggested  by  Liepmann.  Following  on  his  suggestion,  it  is  interesting  to  compare 
that  to  the  Taylor  scale  for  a  turb\ilent  jet,  for  example.  For  turbulence  in  the 
far  field  of  a  jet,  the  Taylor  scale.  At,  can  be  estimated  from  the  Taylor  Reynolds 
number  on  the  jet  axis  (Eq.  12).  This  is  approximately  given  by 

Rct  1.4Re*/2  ,  (17) 

on  the  jet  axis.^°’^^  Using  the  value  of  u'  ~  0.25  Uc,  on  the  axis  of  the  turbulent  jet, 
and  S(x)  ~  0.4  (i  —  lo)  for  the  local  jet  diameter,  we  obtain 

^  ~  2.3Re“'/2  ,  (18) 

which  is  a  little  smziller  but  close  to  the  Liepmann  laminar-layer  thickness,  Al  (a 
prefactor  of  2.3  for  At,  vs.  5.0,  for  Al),  especially  considering  that  it  is  estimated 
from  flow  properties  on  the  jet  axis. 

An  appropriate  inner  viscous  scale,  A„,  can  be  estimated  in  terms  of  the  wave- 
number  k^,,  where  the  energy  spectrum  deviates  from  the  -S/s  power-law  behavior, 
or,  k„XK  —  1/8.®^’®^  This  yields,®^ 


50  Ak  . 


(19) 
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To  estimate  the  Reynolds  ntunber  and  outer  scale  dependence  of  A„,  we  can  use  the 
expression  from  FViehe  tt  for  the  energy-dissipation  rate  on  the  jet  axis,  in 

the  far  field,  t.e.. 


where  Uq  is  the  jet  nozzle  velocity,  da  is  the  jet  nozzle  diameter,  and  xq  is  the  virtual 
origin  of  the  far  field  turbulent  flow.  Substituting  in  Eq.  13  we  then  have 

^  -  0.95  Re-’/^  ,  (21) 

0 

and,  therefore,  for  a  turbulent  jet, 

^  ~  50  Re"®/^  .  (22) 

0 

Substituting  for  Al,  A^,  and  Ak  in  Eq.  16,  we  obtain 

Re-’/*  <  ^  a  50Re-’/*  «  4  <  ^  “  5.0Re-‘/’  <  1  .  (23) 

0  0  0 

The  range  of  intermediate  inviscid  scales,  t.e.,  scales  smaller  than  Al  but  larger 
than  A„,  can  be  seen  to  grow  rather  slowly  with  Reynolds  ntunber.  Specifically,  the 
ratio 

A/  =  4=^  ,  (24a) 

^min 

which  measures  the  extent  of  the  uncoupled  range  of  spatial  scales,  t.e.,  the  number 
of  viscotis  scales  within  a  Taylor  scale,  is  given  by 

M  «  0.1  Re^/*  ,  (24b) 

where  the  (approximate)  factor  of  0.1  was  estimated  for  a  turbulent  jet.  This 
is  indicated  schematically  in  Fig.  11.  In  other  flows,  we  can  expect  the  imcoupled 
range  of  scales  to  exhibit  the  same  Reynolds  number  dependence  (Ek^s.  24a,b),  with, 
possibly,  a  different  prefactor,  however. 


On  the  basis  of  these  observations,  it  can  be  argued  that  a  necessary  condition 
for  fully-developed  turbulence  and  the  1941  Kolmogorov  similarity  ideas  to  apply 
is  the  existence  of  a  range  of  scales  that  are  imcoupled  from  the  large  scales,  on  the 
one  hand,  and  are  free  from  the  effects  of  viscosity,  on  the  other.  Considering  that 


we  must  have 


■^L  ^ 

A(»  Amin 


Af  >  I  , 


(24c) 
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log(R<) 

Fig.  11  Reynolds  number  dependence  of  spatial  scales  for  a  turbulent  jet. 

with  some  margin,  we  see  that  the  existence  of  such  a  range  of  scales  requires  a 
minimiun  Reynolds  number  of  the  order  of  10^  (Eq.  24b);  a  value  in  accord  with  the 
minimum  Reynolds  munber  identified  for  transition  to  fully-developed,  well-mixed 
tmbulent  fiows. 

Jimenez  et  al.  performed  measurements  of  velocity  fiuctuations  in  a  two-dimen¬ 
sional  shear  layer  and  found  a  power-law  re^me  in  the  energy  spectrum,  with  an 
exponent  close  to  -Va,  developing  in  the  neighborhood  of  the  mixing  transition.®® 
Subsequent  investigations  of  the  mixing  transition  by  Huang  and  Ho  also  associated 
the  development  of  a  -  Va  spectral  regime  with  the  mixing  transition,  correlating  it, 
however,  with  the  number  of  pairings  rather  than  with  local  values  of  the  Reynolds 
number.®®  Nevertheless,  in  both  these  investigations,  the  Reynolds  number  in  the 
vicinity  of  the  mixing  transition  and  the  development  of  the  -Va  spectrum  regime 
was  fotmd  to  be  in  the  range  of  3  x  10®  <  Re(x)  <  10®,  in  accord  with  the  range 
documoited  in  Fig.  2.  To  the  extent  that  the  appearance  of  a  -  Va  spectral  regime 
marks  the  development  of  an  inertial  range  of  scales  and  the  applicability  of  the 
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1941  Kolmogorov  ideas,  these  experiments  lend  further  credence  to  the  ansatz. 


6.  Conclusions 

The  preceding  discussion  of  the  experimental  evidence  and  the  theoretical 
ansatz  supports  the  notion  that  fully-developed  turbulence  requires  a  TniniTniiTn 
Reynolds  niimber  of  the  order  of  10'*  to  be  sustained.  This  value  must  be  viewed  as 
a  necessary,  but  not  sufficient,  condition  for  the  flow  to  be  fully-developed.  Presently 
available  evidence  suggests  that  both  the  fact  that  the  phenomenon  occmrs  and  the 
range  of  values  of  the  Reynolds  number  where  it  occurs  are  universal,  ».c.,  indepen¬ 
dent  of  the  flow  geometry. 

On  the  other  hand,  how  sharp  this  transition  is  dots  appear  to  depend  on  the 
details  of  the  flow.  In  particular,  it  is  remarkably  sharp,  as  a  function  of  Reynolds 
number,  in  the  (Couette-Taylor)  flow  between  concentric  rotating  cylinders.  It  is  less 
well-deflned  for  a  shear  layer  and,  among  the  flows  considered,  the  least  well-deflned 
for  turbulent  jets.  Perhaps  an  explanation  for  this  variation  lies  in  the  definition 
of  the  Reynolds  number  itself  (Eq.  1)  and  the  manner  in  which  the  various  factors 
that  enter  are  specified  for  each  flow.  In  the  case  of  the  Couette-Taylor  flow,  for 
example,  both  the  velocity  t/cT  =  H  a  and  the  spatial  scale  ^cT  =  b  —  a,  where  Q 
is  the  difierential  rotation  rate,  with  a  and  b  the  inner  and  outer  cylinder  radii,  are 
well-defined  by  the  flow-boundary  conditions.®* 

In  the  case  of  a  zero  streamwise  pressure-gradient  shear  layer,  the  velocity 
Uai  =  AU  =  Ui  —  U2  is  &  constant,  reasonably  well  specified  by  the  flow  bound¬ 
ary  conditions  at  a  particular  station.  The  length  scale  =  {6{x,t))^, 

however,  must  be  regarded  as  a  stochastic  variable  in  a  given  flow  with  a  relatively 
large  variance.  The  Reynolds  ntunber  for  the  shear  layer  is  then  the  product  of  a 
well-defined  variable  and  a  less  well-defined,  stochastic  variable. 

In  the  case  of  a  turbvilent  jet,  both  the  local  velocity  =  Uj{x)  =  { Uc(a:,  t) 
and  the  length  scale  6j  =  ^j(x)  =  (^j(x,<)),,  or  6j(x)  =  {R{x,t))^,  the  local 
jet  radius,  must  be  regarded  as  stochastic  flow  variables,  each  with  its  own  large 
variance.  The  Reynolds  niunber  for  the  jet  is  then  the  product  of  two  stochastic 
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variables  and,  as  a  consequence,  its  local,  instantaneous  value  is  the  least  well- 
defined  of  the  three. 

Viewing  the  Reynolds  ntunber  itself  as  a  stochastic  variable,  it  would  appear 
that  the  hierarchy  of  the  sharpness  of  the  transition  to  the  fvilly-developed  turbu¬ 
lent  state  is  correlated  with  the  sharpness  with  which  the  fiow  and  the  boundary 
conditions  allow  the  values  of  the  local  Reynolds  number  to  be  specified  to  the 
dynamics. 

A  related  issue  also  arises  as  a  consequence  of  the  definition  of  the  local 
Reynolds  number.  As  noted  in  the  discussion  of  Eqs.  1  and  2,  the  Reynolds  number 
for  a  shear  layer  increases  with  the  downstrezun  coordinate,  whereas  the  Reynolds 
number  for  a  jet  is  a  constant  of  the  fiow.  As  a  consequence,  a  shear  layer  may 
possess  regions  with  local  Reynolds  numbers  below  the  minimum  and  transition  to 
fully-developed  turbulence,  within  the  spatial  extent  of  the  same  flow,  if  its  strezun- 
wise  extent  is  large  enough.  A  turbulent  jet,  on  the  other  hand,  is  either  fully 
developed  over  its  whole  extent,  or  is  not.  This  is  also  relevant  to  the  description 
and  dynamics  in  other  flows. 

As  regards  fully-developed  turbulent  flow,  the  presently  available  evidence  does 
not  support  the  notion  of  Reynolds-nxunber-independent  mixing  dynamics,  at  le2ist 
in  the  case  of  gas-phase  shear  layers  for  which  the  investigations  span  a  large  enough 
range.  In  the  case  of  gas-phase  ttirbulent  jets,  presently  available  evidence  ad¬ 
mits  a  flame  length  stoichiometric  coefiicient  A  (c/.  Ek^.  3)  tending  to  a  Reynolds- 
ntunber-independent  behavior  (c/.  Fig.  7).  We  appreciate,  however,  that  the  range 
of  Reynolds  munbers  spanned  by  experiments  to  date  may  not  be  large  enough  to 
provide  us  with  a  definitive  statement,  at  least  as  evidenced  by  the  range  required 
in  the  case  of  shear  layers  (c/.  Figs.  8).  Secondly,  the  flame-length  virtual  origin,  B 
(Eq.  3),  possesses  a  maximum  in  the  neighborhood  of  the  transition  Reynolds  num¬ 
ber  of  Re  w  2  X  10^  and  does  not  appear  to  attain  a  Reynolds-number-independent 
behavior  in  the  same  range  of  Reynolds  numbers.^^  We  should  also  recall  that 
the  torque  in  Couette-Taylor  flow  does  not  attain  a  Reynolds-ntunber-independent 
behavior  to  the  highest  values  of  the  Reynolds  number  investigated.^®  Neither,  of 
course,  does  the  skin-friction  coefficient  in  a  turbulent  boundary  layer  over  a  smooth 
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flat  plate. 

lu  comparing  shear-layer  with  turbulent-jet  mixing  behavior,  the  more  impor¬ 
tant  conclusion  may  be  that  they  appear  to  respond  in  the  opposite  way  to  Schmidt 
number  effects,  t.e.,  gas-  vs.  liquid-phjise  behavior.  SpecificsJly,  it  is  high-Schmidt 
number  (hquid-phase)  shezu*  layers  that  exhibit  a  low  Reynolds-number  dependence 
in  chemical  product  formation,  if  any  (cf.  Fig.  8b).  In  contrzist,  it  is  gas-phase  turbu¬ 
lent  jets  that  exhibit  an  almost  Reynolds-number-independent  normalized  variance 
of  the  jet-fluid  concentration  on  the  jet  axis,  with  a  strong  Reynolds-number  de¬ 
pendence  foimd  in  hquid-phase  jets,  in  the  same  Reynolds-unber  range  (cf.  Fig.  4). 

To  summarize,  recent  data  on  turbulent  mixing,  as  well  as  evidence  geimered 
in  other  contexts,  support  the  notion  that  fully-developed  turbulent  flow  requires 
a  minimum  Reynolds  number  of  10'*,  or  a  Taylor  Reynolds  number  of  Rct  «  10^, 
to  be  sustained.  Conversely,  turbulent  flow  below  this  Reynolds  ntimber  cannot  be 
regarded  as  fully-developed  and  can  be  expected  to  be  quahtatively  different. 

The  manifestation  of  the  transition  to  this  state  may  depend  on  the  particular 
flow  geometry,  e.g.,  the  appearance  of  streamwise  vortices  and  three-dimensionahty 
in  shear  layers.  Nevertheless,  the  fact  that  such  a  transition  occmrs,  as  well  as 
the  approximate  Reynolds  number  where  it  is  expected,  appears  to  be  a  universed 
property  of  turbulence.  It  is  observed  in  a  wide  variety  of  flows  and  turbulent  flow 
phenomena. 

In  contrast,  studies  of  mixing  in  fully-developed  turbulent  jets  and  shear  layers 
suggest  that  we  cannot  hope  for  a  imiverseil  description  of  turbulent  mixing.  The 
dimensionless  parameters  that  scale  the  relative  importance  of  the  molecular  diffu- 
sivity  coefficients,  such  as  viscosity  and  species  diffusivity,  must  not  only  enter  in 
this  description,  but  are  Ukely  to  do  so  in  a  non-universed  way. 

It  is  interesting  that  transition  to  turbulence,  at  intermediate  values  of  the 
Reynolds  number,  appears  to  be  imiversal,  whereas  mixing  in  fully-developed  tur¬ 
bulence,  at  high  Reynolds  numbers,  does  not. 
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A  method  for  computing  one-dimensional  unsteady  compressible 
flows,  with  and  without  chemical  reactions  is  presented.  This  work  has 
focused  on  the  accurate  computation  of  the  discontinuous  waves  that 
arise  in  such  flows.  The  main  feature  of  the  method  is  the  use  of  an 
adaptive  Lagrangian  grid.  This  allows  the  computation  of  discon¬ 
tinuous  waves  and  their  interactions  with  the  accuracy  of  front-tracking 
algorithms.  This  is  done  without  the  use  of  additional  grid  points 
representing  shocks,  in  contrast  to  conventional  front-tracking 
schemes.  The  Lagrangian  character  of  the  present  scheme  also  allows 
contact  discontinuities  to  be  captured  easily.  The  algorithm  avoids 
interpolation  across  discontinuities  in  a  natural  and  efficient  way.  The 
method  has  been  used  on  a  variety  of  reacting  and  non-reacting  flows 
in  order  to  test  its  ability  to  compute  accurately  and  in  a  robust  way 
complicated  wave  interactions.  S)  i  sss  Academic  Pieu.  inc 


1.  INTRODUCTION 

Several  methods  for  computing  unsteady  inviscid  com¬ 
pressible  flows  have  appeared  in  the  literature  in  recent 
years.  The  emphasis  has  been  on  the  ability  of  these  numeri¬ 
cal  schemes  to  compute  accurately  discontinuous  waves 
which  develop  and  their  interactions. 

High-resolution  shock-capturing  methods  for  hyperbolic 
conservation  laws  is  one  category  of  such  methods  which 
have  been  used  successfully  in  recent  years.  A  basic  feature 
of  these  methods  is  that  the  conservative  formulation  is  used 
which  allows  for  shocks  and  their  interactions  to  be 
captured  automatically  without  special  effort.  This  is 
characteristic  of  all  older  shock-capturing  methods, 
such  as  the  Lax-WendrofT  scheme  [8],  the  MacCormack 
scheme  [  10],  the  original  Godunov  scheme  [5],  In  all  such 
methods,  discontinuous  waves  of  the  solution  are  represen¬ 
ted  as  steep  fronts,  i.e.,  smeared  over  a  finite  number  of  com¬ 
putational  cells.  A  second  and  more  important  feature  of 
recent  high-resolution  schemes  is  the  special  effort  which  is 
made  to  achieve  higher  order  spatial  and  temporal  accuracy 
so  as  to  represent  discontinuous  waves  of  the  solution  as 
accurately  as  possible,  i.e.,  to  reduce  the  smearing  effect 


which  is  typical  of  all  shock-capturing  methods.  Such 
schemes  are  the  TVD  schemes  [6,  7],  the  various  MUSCL- 
type  schemes  [15],  the  PPM  scheme  [3]  (piecewise 
parabolic  method),  etc.  A  comparative  study  of  some  of 
these  schemes  for  real  gases  is  given  in  a  review  article  by 
Montagnc  et  al.  [11].  The  basic  high-resolution  shock¬ 
capturing  methods  have  been  developed  for  nonlinear  scalar 
hyperbolic  conservation  laws.  It  is  for  this  case  that  there 
exists  a  sound  mathematical  theory.  For  nonlinear  hyper¬ 
bolic  systems  of  equations  in  one  space  variable  the  theory 
is  not  as  clear  and  the  numerical  methods  used  for  these 
systems  apply  formally  the  same  techniques  as  in  the  scalar 
case,  but  with  the  additional  use  of  exact  or  approximate 
Riemann  solvers.  A  classical  Riemann  problem  is  solved 
locally  at  each  computational  cell  boundary  in  order  to 
compute  the  various  flux  terms  required.  This  is  the  essen¬ 
tial  ingredient  of  the  original  Godunov  scheme  and  it  is  pre¬ 
sent  in  most  successful  high-resolution  schemes.  The  various 
flux-vector  splitting  techniques  [14,  16]  have  essentially 
incorporated  in  them  an  approximate  Riemann  solver. 
Finally,  their  extension  to  more  than  one  space  dimension 
is  usually  done  by  treating  each  spatial  dimension 
separately. 

Another  category  of  numerical  schemes  that  have  been 
used  is  that  of  the  shock-fitting  or  front-tracking  methods. 
Although  they  have  not  been  used  as  extensively  as  the 
shock-capturing  methods,  they  have  been  quite  successful  in 
one-dimensional  problems.  A  good  review  of  these  methods, 
as  well  as  of  many  shock  capturing  methods,  is  given  by 
Moretti  [12].  These  schemes  are  typically  based  on  a  non¬ 
conservative  formulation  and  an  effort  is  made  to  detect  and 
identify  the  various  discontinuous  waves  and  compute  their 
interactions  explicitly.  This  is  usually  accomplished  by 
introducing  additional  computational  elements  repre¬ 
senting  such  waves  and  using  the  Rankine-Hugoniot  jump 
conditions.  This  technique  leads  to  complex  programming 
logic.  Identifying  the  waves  and  computing  their  interac¬ 
tions  accurately  is  crucial  for  obtaining  a  meaningful  and 
stable  solution.  For  flows  with  complicated  wave  inter- 
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actions  such  schemes  may  not  be  as  robust  as  the  shock¬ 
capturing  schemes,  even  in  one  space  dimension. 

The  research  presented  in  this  paper  is  part  of  a  greater 
effort  which  aims  to  combine  the  characteristics  of  the 
above  two  categories  of  numerical  schemes  and  to  develop 
a  method  which  will  share  the  advantages  and  eliminate 
most  of  the  disadvantages  of  both.  This  has  been  accom¬ 
plished  in  the  case  of  one-dimensional  flow  by  the  scheme 
presented  in  this  paper.  The  increased  accuracy  which  is 
provided  in  the  computation  of  complicated  wave  inter¬ 
actions  and  its  robustness  have  made  this  scheme  especially 
valuable  for  the  computation  of  reacting  gas  flows,  where 
detonation  waves  are  present. 

The  scheme  is  based  on  a  conservative  shock -capturing 
Godunov-type  scheme,  very  much  like  van  Leer’s  MUSCL 
scheme  [15].  The  new  feature,  introduced  here,  is  an  adap¬ 
tive  Lagrangian  grid  which  increases  the  accuracy  with 
which  discontinuous  waves  and  their  interactions  are 
computed.  Without  introducing  additional  computational 
elements,  i.e.,  refining  the  grid,  or  special  computational 
elements  to  represent  these  waves,  the  shocks  and  contact 
discontinuities  are  computed  as  true  discontinuities, 
without  the  smearing  effect  typical  of  shock-capturing 
methods.  This  makes  the  scheme  different  from  adaptive 
mesh  refinement  schemes  (e.g.,  see  Berger  and  Oliger  [2]), 
which  smear  discontinuities,  although  on  a  much  finer  local 
grid.  The  basic  conservative  shock -capturing  capabilities  of 
the  scheme  are  not  diminished.  The  scheme  is  endowed  with 
the  capability  to  track  various  fronts  and,  thus,  the  shock¬ 
capturing  and  the  front-tracking  ideas  are  combined 
properly.  It  is  important  to  note  that  the  adaptive  grid 
strategy,  to  a  certain  degree,  is  independent  of  the  particular 
solver.  Any  Godunov-type  scheme  may  be  used.  The 
Riemann  solver  is  the  link  that  provides  the  information 
about  local  wave  interactions  needed  for  the  adaptive 
procedure. 

It  was  deemed  interesting  to  try  this  scheme  on  one¬ 
dimensional  flows  of  reacting  gases  in  light  of  the  increased 
accuracy  and  robustness  with  which  detonation  waves  and 
their  interactions  could  be  computed.  The  interest  in  such 
flows  is  evident  by  the  number  of  papers  appearing  in 
the  literature.  For  example,  numerical  calculations,  with 
increased  accuracy,  of  the  one-dimensional  instability  of 
plane  detonation  waves  may  be  of  great  interest  in  confirm¬ 
ing  existing  theories  which  are  based  on  linear  stability 
analysis  (e.g.,  see  Lee  and  Stewart  [9]).  The  present  scheme 
is  able  to  reduce  the  error  caused  by  the  numerical  smearing 
of  the  leading  shock  of  the  detonation  wave.  This  error 
may  be  very  important  in  the  development  of  detonation 
instability. 

The  computer  code  developed  is  also  able  to  compute 
one-dimensional  cylindrically  and  spherically  symmetric 
flows,  as  well  as  plane  flows  with  area  change.  It  is  thus 
possible  to  compute  explosions  and  implosions  and  study 


the  effect  of  curvature  on  detonation  wave  speed  and 
stability.  Most  of  the  results  presented  are  basically  valida¬ 
tion  runs  and  calculations  demonstrating  the  abilities  of  the 
method  and  the  potential  use  for  specific  one-dimensional 
problems  of  interest.  All  results  shown  are  for  a  perfect  gas. 
The  difliculty  of  incorporating  a  general  equation  of  state  is 
the  same  as  in  most  schemes  and  independent  of  the  main 
feature  of  the  present  scheme,  i.e.,  the  adaptive  Lagrangian 
grid  strategy. 


2.  NUMERICAL  METHOD 

2.1.  Mathematical  Formulation 

The  inviscid  flow  of  a  reacting  mixture  of  calorically 
perfect  gases  is  considered.  The  assumption  of  a  simplified 
reacting  mixture  is  made,  according  to  which  there  are  two 
species  present  at  any  time,  the  reactant  and  the  product. 
The  reactant  is  converted  to  the  product  by  a  one-step  irre¬ 
versible  exothermic  chemical  reaction.  This  assumption  is 
made  in  order  to  compare  with  the  many  theoretical  and 
numerical  results  which  are  available  in  the  literature  for 
this  case.  The  chemical  reaction  rate  is  given  by  the  standard 
Arrhenius  law 


i= -.Kzr'expf -£//?,  D,  (la) 


where  z  is  the  mass  fraction  of  unburnt  gas,  £  is  a  positive 
constant,  which  essentially  gives  a  time  scale,  E  is  the  activa¬ 
tion  energy  of  the  chemical  reaction,  Rg  is  the  gas  constant, 
T  is  the  absolute  temperature,  and  a  is  also  a  constant.  The 
simplified  Arrhenius  model,  where  the  reaction  rate  is  a  step 
function  depending  on  the  temperature,  has  al.''>  been  used. 
For  the  simplified  model  the  rate  is  given  by 

z= -KzH(T-T,),  (lb) 


where 


H(x) 


1,  x>0, 

0,  X  ^  0, 


(2) 


and  Tc  is  a  given  critical  temperature. 

The  problem  under  consideration  is  a  special  case  of  the 
general  problem  of  solving  numerically  the  nonlinear 
hyperbolic  system  of  the  form 


dF(U) 
dt  ^  dx 


=Gm 


(3) 


where  1/  is  the  appropriate  solution  vector.  As  usual,  x 
denotes  the  Eulerian  space  variable.  If  the  Lagrangian 
formulation  is  used,  a  system  of  exactly  the  same  form  is 
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obtained.  The  space  variable  jc,  then,  is  replaced  by  a 
Lagrangian  space  variable  and  the  flux  vector  F\U)  is 
changed  appropriately.  For  non-reacting  flow,  G(U)  =  0. 
Most  numerical  methods  use  Eq.  (3)  as  their  starting  point 
and,  using  a  finite  volume  discretization,  obtain  the  scheme 
of  the  following  general  form 

./2  -  ./2) + Gy.  (4) 

giving  the  solution,  in  an  average  sense,  in  the  y  th  cell  at  the 
time  level  n  -i- 1.  The  numerical  flux  terms  Fare  computed  at 
the  boundaries  of  each  cell.  An  important  feature  of  every 
numerical  method  is  the  calculation  of  these  flux  terms  in  a 
way  that  guarantees  stability  and  high-order  accuracy. 

A  slightly  different  approach  will  be  taken  in  deriving  the 
present  scheme.  Eventually,  it  will  be  of  the  general  form 
given  in  Eq.  (4).  It  is  useful  to  formulate  the  problem  by 
writing  the  conservation  laws  in  integral  form  for  an 
arbitrary  control  volume  F(/),  whose  bounding  surface  S(/) 
moves  with  a  velocity  Ui,  (Reynolds’  transport  theorem). 
These  equations  will  be  applied  to  each  computational 
volume  of  the  discrete  numerical  scheme.  This  is  done  so 
that  the  conservation  equations  and  their  discrete  counter¬ 
parts  are  written  in  a  way  which  is  independent  of  the 
Eulerian  or  Lagrangian  formulation  that  will  be  adopted 
eventually.  Moreover,  it  is  easier  to  see  from  these  equations 
how  the  idea  for  the  adaptive  nature  of  the  grid  is  motivated. 
The  conservation  equations  in  integral  form  are 


equations,  e^  is  the  total  specific  energy,  which  includes  the 
chemical  energy,  i.e., 

e^  =  e  +  ku^  +  qo2,  (9) 

where  e  is  the  specific  internal  energy,  qo  is  the  heat  release 
of  the  chemical  reaction,  and  u  =  |u|  is  the  magnitude  of  the 
fluid  velocity.  The  perfect  gas  assumption  is  also  made,  i.e., 

p  =  {y-^)pe.  ( 10) 

Since  the  boundaries  of  the  computational  cells  will  be 
moving,  it  is  useful  to  consider  the  flow  map 

x  =  X(4,0,  (11) 

which  gives  the  position  of  the  fluid  particle  that  was 
initially  (/  =  0)  at  the  position  Thus,  4  is  convenient 
Lagrangian  marker  for  the  fluid  particles  in  the  flow.  If  the 
Lagrangian  approach  is  taken,  the  local  boundary  velocity 
is  equal  to  the  local  fluid  velocity,  i.e.,  U),  =  u  in  Eqs.  (5)-(8). 

2.2.  Spatial  and  Temporal  Discretization 

Consider  now  the  case  of  one-dimensional  flow.  A  finite 
volume  formulation  is  used,  i.e.,  space  is  discretized  by  a  set 
of  computational  cells  as  shown  in  Fig.  1.  The  conservation 
equations  are  now  written  for  the  yth  cell  of  the  computa¬ 
tional  grid 


[  pdVJrl  p(u-Ub)</S  =  0,  (5) 

dt  Jy(r)  ■'S(i) 

d 


[  pudy+ 

pu(u-Ub)dS 

'ni) 

S(l) 

+  f  pdS 

=  0, 

(6) 

•'S(l) 

f  pe,dy+ 

[  pc,(u  -  Ub)  •  </S 

V(/) 

■'siri 

-1-  [  pu  dS  =  0, 

(7) 

■'S(l) 

[  pz  dV  +  I 

pr(u-Ub)dS 

S(0 

-f  zpdV=0. 

(8) 

dm: 

+  (P  dWb)y+  1/2  -  (P  dMb)>-  1/2  =  0. 


—  (mjUj)  +  {pu  dMb)>+  1/2  -  (P“  dUb)>-  1/2 

+  Py  +  1/2  “  Pj-  1/2  = 
d 

—  (m^e^j)  (pc,  d«b)y+ 1/2  -  (P^I  dUb);_  1/2 

+  (“P  )y  + 1/2 -{«p)y- 1/2  =  0, 
d 

^("'y‘>)+(PZdMb)y+i/2 
-(prdMb)y-i/2-z>"»y  =  0, 


j*cell 


space 

variable 


(12) 


(13) 


(14) 


(15) 


FIG.  I.  Finite  volume  discretization  in  one  space  dimension.  The 
These  are  written  for  an  arbitrary  control  volume  F(l),  space  variable  can  be  the  Eulerian  x  or  the  Lagrangian  The  boundaries 
whose  bounding  surface  S(t)  has  a  velocity  Ub.  In  the  above  of  theyth  cell  are  denoted  by  the  subscriptsyi ). 
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where 


Ill 

a 

«-«b. 

mjm 

(••»/+ 1/2 

1  pdx. 

r-'y+  1/2 

pu  dx. 

1  P(^+; 

1/2 

r  't,  *  1/2 

nt: 


'l  I, -2 


(16) 


Average  values  of  all  quantities  in  the yth  cell  are  denoted 
by  the  subscript  j  and  values  of  various  quantities  at  the  two 
boundaries  of  the  cell  are  denoted  by  the  subscripts  j±{. 
Note  that  average  values  are  mass-averaged  values.  By 
defining 

F„  =  pAuy„ 

=  pu  dUb  +  P. 

(17) 

F,  s  pe,  dUb  -t-  pu, 

F.spzAui,, 


the  equations  of  motion  can  be  written  in  the  more  familiar 
form 


The  equation  of  state  (10)  provides  the  means  for  com¬ 
puting  the  average  pressure  in  the yth  cell, 

=  (21) 

The  exact  integral  conservation  laws  have  been  written 
for  each  computational  cell.  Equations  (18)  will  now  be 
integrated  in  time  explicitly.  The  basis  for  the  method  is  a 
conservative  Godunov-type  scheme  similar  to  the  MUSCL 
scheme  introduced  by  van  Leer  [15].  The  procedure 
followed  in  solving  these  equations  is  similar  to  that  used  in 
most  methods,  which  are  higher-order  extensions  of  the 
original  Godunov  scheme.  At  every  time  instant,  average 
values  of  the  solution  are  known  in  each  computational  cell. 
Linear  variations  of  the  primitive  variables,  i.e.,  density  p, 
pressure  p,  and  velocity  u,  are  assumed  in  each  cell.  A 
Riemann  problem  is  then  set  up  locally  at  each  cell  interface. 
The  solution  to  this  problem  gives  the  velocity,  pressure  and 
density  needed  to  compute  the  flux  terms  (17).  The  different 
feature  in  the  present  scheme  is  that  the  Lagrangian 
formulation  is  used  instead  of  the  Eulerian  and  that  an 
adaptive  grid  is  used. 

So  far,  the  fact  that  the  Lagrangian  formulation  is  being 
used,  has  not  appeared  explicitly  in  the  description  of  the 
method.  It  is  now  that  this  choice  is  made  and  all  quantities 
are  considered  as  functions  of  time  t  and  the  Lagrangian 
space  coordinate  The  interpolation  procedure  is  carried 
out  in  i^-space  and,  assuming  linear  variation,  the  generic 
quantity  q  varies  as 


-^+(fJy../2-(7^™)y-./2  =  0, 

d 

—  (nijUj)  +  (F„)j^  1/2  -  (F„)j_  1/2  =  0, 
d 

—  (mje,j)  +  (F,)j^  ,/2  -  (F,)y- 1/2  =  0, 

d 

—  (nijZj)  +  (F,)j^  ,/2  -  (F,)j_  ,/2  =  nijZj 


(18) 


Note  the  extra  degree  of  freedom  provided  in  the  flux  terms 
by  the,  as  of  yet  unspecified,  term  Au^,.  The  motion  of  the 
cell  boundaries  is  determined  by 


Chilli -if.  \ 

-l“b;y±i/2 


and  the  average  density  in  each  cell  is  given  by 

Pi^^' 


(19) 


m. 


’fy+  1/2  ~^j  -  1/2 


(20) 


9(0  =  9y+(9i)y(^“<?y),  (22) 

in  the yth  cell,  where  qj  is  the  mass-averaged  value  in  the  cell, 
is  the  center  of  the  cell  (in  Lagrangian  space),  and  (q^)j  is 
the  slope  of  q  in  this  cell,  which  is  assumed  to  be  constant. 
Note  that  discontinuities  of  these  quantities  are  allowed  at 
the  cell  interfaces,  as  shown  in  Fig.  2. 

The  slopes  are  chosen  using  the  van  Leer  slope  limiter 
[17],  but  the  adaptive  nature  of  the  grid,  which  will  be 


q  A 
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‘‘j 

n 

9/1 

_ 1 

j-3/2 

j-  1/2 

jiKrell 

j+l/2 

Space 
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FIG.  2.  Linear  variation  of  the  generic  quantity  q  in  the  y'th  cell. 
In  general,  q  is  discontinuous  at  the  cell  interfaces. 
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described  next,  makes  the  choice  of  limiter  less  important 
than  in  the  typical  higher-order  Godunov-type  schemes.  In 
fact,  the  adaptive  grid  allows  more  freedom  in  choosing  the 
interpolation  scheme,  because  additional  information  on 
the  location  of  the  various  waves  is  always  available  at 
each  time  instant.  The  slope  (q^)^  is  computed,  following 
van  Leer  [17],  by 


(qi)^  =  &ve{q^,q^). 


where 


<]i  = 


ave(.v,  >’)s 


.V  +  >’ 
2 


(-V-V)-  I 
x^-  +  y^  +  c^y 


(23) 


(24) 

(25) 


and  c-  is  a  small  constant  (c^  <1 1 ). 

At  each  cell  interface,  two  constant  states  q  and  q  *  are 
required  to  be  used  as  the  initial  condition  for  the  Riemann 
problem.  There  are  many  ways  of  doing  this.  One  way  is  to 
specify  for  the  7  +  5  interface 


where  po  is  the  initial  density  (/  =  0)  and  a  is  the  speed  of 
sound.  The  constant  states  q^  are  then  determined  as  the 
averages  of  Eq.  (22)  over  the  domains  Fig.  3.  These 
domains  correspond  to  the  full  timestep  di.  This  is  equiv¬ 
alent  to  tracing  the  characteristics  back  from  the  time 
t-yAtjl  and  using  the  linear  profile  (22).  This  ensures 
second-order  accuracy  in  time. 

The  discrete  scheme,  giving  the  solution  at  the  time  level 
n  -F  1  from  the  solution  at  the  previous  time  level  n,  can  now 
be  written  as 

,,,  -  {FJj_ 

=  (OT^Cy)"-  J/[(F,),+  1/2-  (/■,);_  1/2],  (28) 

(m,-)’’"  '  =  {m,=y  -  ^/[(  A),,  ,  2  -  (A),^  ,,2] 

^J±  M2  ~  ^i±  1,2  "I"  '4/(«|,)y+  1/2, 


9/ +1/2  4"  (^ <);  (<Jy+ 1/2  ^y  )i 

/ 

9y*+  1,'2  “  ?/  +  I  +  (94)/+  1  (<Jy+  1/2  ~  <?/+  1  ). 

i.e.,  the  values  of  q  on  cither  side  of  the  interface,  as  given  by 
Eq.  (22).  Using  these  states  does  not  ensure  second-order 
accuracy  in  time.  The  method  used  in  the  present  scheme  is 
shown  in  Fig.  3.  The  domain  of  dependence  of  <J  =  (Sy+i/2 
over  the  time  interval  At  is  estimated  by  the  characteristics 
at  the  time  level  t.  In  the  Lagrangian  formulation  of  the 
problem  the  characteristic  speeds  are  given  by 

c±  =  ±-^a,  (27) 

Po 


j+l/2 

FIG.  3.  The  constant  states  q*,  which  are  to  be  used  as  the  initial 
condition  for  the  Riemann  problem  at  the  interface  j+\,  are  obtained  by 
averaging  the  linear  interpolant  over  the  domains  of  dependence  Ji*. 
These  domains  correspond  to  the  full  timestep  d(. 


where  the  numerical  fluxes  F„,  F^.  F^,  and  A  are  given  by 
Eqs.  (17),  using  the  solution  of  the  Riemann  problem.  The 
average  boundary  velocity  u^,  for  each  interface  is  still 
unspecified,  but  for  the  majority  of  interfaces  «b  =  «  and  the 
last  of  Eqs.  (28)  is  second-order  accurate  in  time.  The 
source  term  in  the  species  equation  is  shown  in  Eqs.  (28)  as 
being  evaluated  at  the  time  level  n.  It  is  better  to  integrate 
the  source  term  in  a  “split”  manner,  i.e.,  integrate  the  first 
four  equations  in  (28)  without  the  source  term  and  use  this 
intermediate  state  to  estimate  the  term  This  splitting 
has  been  implemented  in  the  present  scheme. 

The  stability  requirement  on  the  timestep  is  that  of  a 
MUSCL  scheme  in  the  Lagrangian  formulation.  No  addi¬ 
tional  stability  problems  arise  due  to  the  adaptive  grid 
strategy  presented  in  the  next  section. 

2.3.  Adaptive  Grid 

The  motivation  for  the  adaptive  grid  comes  from  the 
definition  of  the  flux  terms,  as  given  by  Eqs.  (17).  The  term 
Jwb,  or,  equivalently,  the  velocity  of  the  cell  boundary  Mj,, 
is  unspecified.  The  idea  is  to  specify  it  at  each  cell  interface, 
so  that  all  important  discontinuous  waves  coincide  with  cell 
boundaries,  at  every  discrete  time  level.  The  solution  of  the 
Riemann  problem  at  a  given  interface  provides  all  the 
information  needed  to  identify  all  the  important  waves 
emanating  from  this  interface,  as  well  as  their  strengths  and 
speeds.  This  information  is  enough  to  specify  Au^,.  Since  all 
important  waves  coincide  with  cell  boundaries,  it  is  guaran¬ 
teed  that,  at  subsequent  time  instants,  the  evolution  of  these 
waves  will  be  determined  properly  by  the  solution  of  the 
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local  Riemann  problems.  In  the  numerical  experiments 
carried  out,  shock  waves  computed  by  the  local  Riemann 
solvers  were  considered  important  enough  to  track  if  the 
shock  Mach  number  was  greater  than  1.01  and  contact  dis¬ 
continuities  were  considered  important  if  the  ratio  of  the 
densities  on  either  side  was  greater  than  l.OS.  These 
parameters  are  quite  conservative.  One  may  want  to  track 
only  the  very  strong  shock  waves  in  the  flow. 

The  grid  is,  basically,  Lagrangian,  i.e.,  most  cell  bound¬ 
aries  move  with  the  local  fluid  velocity  and,  hence  Juf,  =  0. 
It  is  easy  to  see  that  the  same  ideas  on  the  adaptivity  of  the 
grid  can  be  used  on  a  grid  that  is  primarily  Eulerian.  The 
same  equations  can  be  applied  directly. 

An  example  of  this  adaptive  procedure  is  shown  in  Fig.  4. 
A  strong  shock  wave  moving  to  the  right  is  computed  by  the 
Riemann  solver  at  the  interface  /  —  5  at  time  /.  The  decision 
is  made  to  assign  a  velocity  to  the  adjacent  cell  boundary 
i+j,  so  that  at  time  t  +  Jt  the  shock  coincides  with  the 
interface  /-(-j.  Another  possibility  would  be  to  have  the 
interface  /  —  |  move  with  the  shock.  The  decision  is  made 
depending  on  which  interface  would  be  required  to  move  a 
shorter  distance  in  Lagrangian  space.  The  shock  speeds  are 
assumed  constant  over  the  time  interval  Jt.  It  is  obvious 
that  the  local  expansion  waves  can  be  tracked  in  the  same 
way.  This  was  not  implemented  in  the  present  scheme, 
simply  to  reduce  the  complexity  of  the  programming. 

It  is  evident  from  this  example  that  a  relation  between  the 
velocities  in  real  space  and  the  velocities  in  Lagrangian 
space  is  needed  to  upxlate  the  Lagrangian  grid.  Consider  the 
motion  of  a  cell  boundary  given  by  the  trajectory  x  = 

This  boundary  is  moving  with  a  velocity  «b  =  ib(t).  which, 
in  general,  is  different  from  that  of  the  fluid  u.  This  motion 
corresprands  to  a  motion  in  ^-space  given  by  the  trajectory 
^  =  ^b(0  with  velocity  Ub  =  ^b(0-  The  relation  between  the 
two  velocities  is  found  with  the  use  of  the  flow  map 

x  =  X(lt),  (29) 


time 


FIG.  4.  The  appropriate  velocity  is  assigned  to  the  cell  interface  i  +  j 
in  order  to  intercept  the  shock  at  the  subsequent  discrete  time  level. 


which  is  essentially  Eq.  (11)  written  here  for  one-dimen¬ 
sional  flow.  The  cell  boundary  motion  is  given  by 


Xb(r)  =  ^(U0,/)  (30) 

and,  hence. 


The  derivative  of  the  flow  map  is  numerically  approximated 
and  assumed  constant  in  each  cell,  i.e.. 


/  ^  ^j+ 1/2  ~  1/2 

C>+ 1/2  ~  1/2 


(32) 


The  velocities  of  the  various  waves,  which  are  computed  by 
the  Riemann  problems,  can  be  translated  into  velocities  in 
(J-space  by  using  Eq.  (31b).  The  Lagrangian  position  of 
each  interface  is  updated  by 


^;±./2=^;±..2+'^/(»’b)>±./2-  (33) 

The  solution  to  the  Riemann  problem  at  each  interface 
provides  sufficient  information  for  the  adaptive  strategy. 
Using  the  exact  Riemann  solver  at  every  interface  is  very 
costly.  To  reduce  the  cost,  various  criteria  were  found  to 
identify  the  cell  interfaces  where  a  strong  discontinuous 
wave  is  suspected  to  be  present,  before  solving  the  Riemann 
problem.  These  interfaces  are  flagged  as  critical  interfaces. 
The  ratio  has  proven  useful  in  detecting  develop¬ 

ing  shocks  in  the  flow.  Where  the  flow  is  smooth,  without 
steep  gradients,  the  above  ratio  is 


s:  1. 


(34) 


The  regions,  where  this  ratio  deviates  from  unity  by  more 
than  10%,  are  considered  critical  regions.  The  full  nonlinear 
Riemann  solver  is  used  only  in  these  regions.  Everywhere 
else  the  simple  acoustic  approximation  to  the  Riemann 
problem  solution  is  used.  It  was  found  in  all  the  numerical 
experiments  performed  that,  in  addition  to  the  above 
criterion,  finding  local  extrema  of  the  slopes  in  pressure, 
density,  and  velocity  was  very  useful  in  determining  these 
regions.  Other  criteria  may  also  be  used.  It  is  important  that 
the  criteria  be  conservative  enough,  so  that  no  critical 
regions  are  missed,  but  they  are  not  crucial  in  detecting 
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discontinuous  waves.  The  detection  of  important  discon¬ 
tinuous  waves  is  ultimately  done  by  examining  the  solutions 
to  the  local  Riemann  problems. 

No  more  than  two  adjacent  critical  interfaces  are  allowed 
at  any  given  time.  In  the  smooth  compression  regions,  the 
interfaces  at  the  maxima  of  the  pressure  gradients  are 
considered  critical.  The  pressure  gradients  are  estimated 
using  simple  finite  differences.  If  the  Riemann  solver  at  these 
locations  computes  sufficiently  strong  discontinuous  waves, 
then  they  are  tracked.  The  critical  interfaces,  which  carry 
these  waves,  are  treated  the  same  way  at  the  next  time  level 
along  with  other  possible  critical  interfaces  that  may  be 
detected.  When  there  are  two  adjacent  critical  interfaces,  the 
two  Riemann  problems  are  solved  simultaneously.  At  this 
point  it  is  decided  if  collisions  will  occur  within  the  timestep 
At.  There  are  six  waves  resulting  from  the  two  Riemann 
problems  and  there  are  four  cell  interfaces  available  to  do 
the  tracking.  The  strongest  waves  are  tracked  and  the  others 
are  ignored.  This  procedure  has  proven  to  be  very  robust  in 
handling  all  possible  wave  interactions. 

Collisions  and  reflections  from  walls  can  be  treated  in  a 
straightforward  way  using  this  adaptive  grid.  A  typical  colli¬ 
sion  case  is  shown  in  Fig.  5.  At  time  t,  two  strong  shock 
waves  at  the  interfaces  i  —  5  and  1  +  5,  are  moving  at  each 
other  with  speeds  that  allow  for  a  collision  before  time 
i-f-d/.  The  Riemann  problems  at  the  interfaces  «  — 5  and 
i  -I-  5  are  solved  at  time  t,  simultaneously.  The  solution 
indicates  that  there  will  be  a  collision  within  the  time  inter¬ 
val  At.  The  time  step  is  adjusted  locally,  i.e.,  only  for  the 
three  cells  /  -  1,  and  /-f- 1,  so  that  at  the  intermediate  time 
instant  the  collision  point  coincides  with  the  cell  boundary 
j-f-j.  The  Riemann  solver  at  this  interface,  at  the  inter¬ 
mediate  time  instant,  will  compute  the  two  shock  waves 
emerging  from  the  collision  and  the  adjacent  cell  boundaries 


FIG.  5.  The  typical  collision  of  two  shocks  is  shown.  The  time  step  is 
adjusted  locally  so  that  the  collision  point  coincides  with  the  cell  interface 
I  +  i  at  the  intennediate  time  step. 


will  be  able  to  track  them  in  the  same  way  at  subsequent 
times.  The  fluxes  at  the  interfaces  /  —  §  and  /  +  5  are  held 
constant  for  the  whole  timestep  At.  This  leads  to  a  robust 
way  of  handling  wave  interactions,  without  loss  of  accuracy. 

3.  RIEMANN  SOLVER 

The  Riemann  solver  is  an  important  ingredient  of  the 
numerical  scheme.  It  provides  the  means  for  computing  the 
velocity  and  the  pressure  at  the  cell  interfaces  and,  thus,  the 
various  flux  terms  required.  It  also  gives  valuable  informa¬ 
tion  about  the  local  waves  emanating  from  each  cell  inter¬ 
face.  As  explained  in  the  previous  section,  the  Lagrangian 
grid  adapts  in  such  a  way  that  important  discontinuous 
waves  and  collision  points  coincide  with  cell  boundaries  at 
each  time  instant.  It  is,  therefore,  necessary  to  be  able  to 
identify  the  waves  emanating  from  these  critical  cell  bound¬ 
aries  at  subsequent  times.  This  is  what  the  Riemann  solver 
accomplishes.  A  variety  of  exact  and  approximate  Riemann 
solvers  have  appeared  in  the  literature  in  recent  years.  In  all 
these  solvers  the  focus  is  on  computing  the  velocity  and 
pressure  of  the  contact  discontinuity,  which  appears  after 
the  breakup  of  the  initial  discontinuity  of  the  Riemann 
problem.  In  the  present  scheme  it  is  crucial  to  identify  the 
exact  wave  pattern  as  well.  This  information  is  used  to 
assign  the  appropriate  velocities  to  adjacent  cell  boundaries 
so  that  all  important  waves  are  tracked  and  to  adjust  the 
time  step  locally  so  that  collisions  are  computed  accurately. 
Moreover,  the  fluxes  at  an  interface  need  to  be  computed 
along  the  ray  =  see  Eq.  {31b).  Most  interfaces  are 
Lagrangian  and  hei.ee,  Vy,  =  0. 

3.1.  Non-reacting  Perfect  Gas 

Consider  the  case  of  the  Riemann  problem  for  inviscid 
flow  of  a  perfect  gas  without  chemical  reactions.  The  initial 
condition  at  time  t  =  0  consists  of  two  constant  states 
denoted  by  the  subscripts  r  and  /.  Note  that  it  is  possible  to 
have  two  different  perfect  gases  on  either  side  of  the  (J  =  0 
location,  as  indicated  by  the  different  specific  heat  ratios, 
i.e.,  fr  and  y,;  see  Fig.  6.  The  space  variable  ^  is  the 


P,.p,.u,,Y, 

Pr'Pr-“r-7r 

k 

FIG.  6.  Initial  condition  for  the  Riemann  problem.  The  variable  i  is 
the  Lagrangian  space  coordinate. 
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Lagrangian  space  coordinate.  At  time  /  =  0*  the  general 
wave  pattern  shown  in  Fig.  7  will  develop. 

There  is  a  wave  moving  to  the  right  (positive  (J)  denoted 
by  R,  a  wave  moving  to  the  left  (negative  O  denoted  by  L, 
and  a  contact  discontinuity  C  which  remains  at  ^  =  0  for  all 
time,  i.e.,  moves  with  the  local  iluid  velocity.  The  waves  R 
and  L  are  either  shocks  or  expansion  waves,  depending  on 
the  initial  condition.  Across  the  contact  discontituity  C  the 
pressure  py  and  velocity  Uf  are  continuous,  but  the  density 
has  a  jump  discontinuity  at  =  0  for  all  time.  The  density  is 
Py,  for  ^  <  0  and  py,  for  ^  >  0.  It  is  known  that  the  solution 
to  this  initial  value  problem  exists  and  is  unique  for 
arbitrary  initial  conditions.  Moreover,  the  solution  is  self¬ 
similar  and  the  shock  waves  propagate  with  a  constant 
velocity  and  strength.  That  is  why  they  are  represented  by 
straight  lines  in  the  ((^,  t)  diagram. 

There  are  four  wave  patterns  possible  for  tins  problem. 
The  solution  will  be  found  for  each  of  these  wave  patterns 
for  the  special  case  of  a  perfect  gas. 

(i)  L-shock,  R-shock.  Across  the  shock  R  the  following 
relationship  holds: 


M,  is  the  shock  Mach  number  defined  by 


Ur-U, 

Or 


(35) 


(36) 


where  I/,  is  the  shock  velocity  and  a,  =  y  y,  PrIPr  is  the 
speed  of  sound  in  the  undisturbed  region  /  .  Similarly  for  the 
shock  L, 


U/-Uy^  2 

a,  y,+  l 


(37) 


where  M,  is  the  shock  Mach  number  for  L.  Equations  (35) 
and  (37)  can  be  solved  for  the  shock  Mach  numbers  to  give 


The  pressure  ratios  across  the  two  shocks  are  given  by  the 
equations 


(40) 

pr  y,+  \ 

(41) 

Pi  Vz+l 

From  Eqs.  (38)-(41 )  one  obtains  a  single  equation  for  the 
unknown  r. 


This  equation  is  solved  numerically  using  a  Newton- 
Raphson  method.  Once  r  is  obtained,  all  other  quantities  of 
interest  follow  from  Eqs.  (35)-(41 ).  The  densities  are  deter¬ 
mined  by 


FIG.  7.  General  wave  pattern  resulting  from  the  breakup  of  the 
original  discontinuity  of  the  Riemann  problem.  C  is  a  contact  discon¬ 
tinuity.  L  and  R  can  be  either  shocks  or  expansion  waves. 


Pp  (y,+  l)Mf 
p,  2-f  (y,-l)A/?’ 


i  =  r,  /. 


(43) 


It  is  important  to  be  able  to  determine  if  this  wave  pattern 
will  develop  for  a  given  initial  condition.  For  this  solution  to 
be  possible,  certain  compatibility  conditions  must  hold. 
These  are  easily  found  by  noticing  that  in  Eqs.  (35)  and  (37) 
the  shock  Mach  numbers  M,  and  M,  must  be  greater 
than  1.  It  then  follows  that  the  following  compatibility 
condition 


U,  <  Uy  <  U, 


(4;) 
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must  hold  or,  equivalently, 


u,  -  u,  ^  0, 


(45) 


(»)  L-expansion  wave,  R-shock.  In  this  case,  Eqs.  (35), 
(36),  (38),  and  (40),  derived  previously,  still  hold  for  the 
shock  wave  R.  In  addition  to  these  equations,  the  following 
equation  gives  the  pressure  ratio  across  the  isentropic 
expansion  wave 


(iv)  L-expansion  wave,  R-expansion  wave.  In  this  case, 
there  are  two  expansion  waves  and  the  pressure  ratios 
across  them  are  given  by 


Pi  ' 

2  a,) 

(50) 

^  =  1 

Pr  ' 

<  2  a,) 

From  Eqs.  (50)  the  following  single  equation  in  r  is 
obtained. 


P- 


1  + 


I  f 


y,-  I  u,  " 

2  a,  } 

y,-lii,-i/,  y,-\a. 

2  u,  2  a') 


II 


(46) 


where  t  n  defined  in  Eq  (35 1  Combining  Eqs.  (38),  (40), 
and  1 46 1,  a  single  equation  in  r  is  obiainnl,  as  in  the 
previous  case. 


Flo 


V.  V'  I  . 

I  •  V - -  f’  * 

4 


P 

P. 


(' 


II 

-0. 


(47) 


which  IS  sciivcd  nucncmaOs  The  dengues  arc  determined 
hv  f  q  1 4 1 1  acrifVM  the  shock  aitd  hy  the  neniropK  relation 


S(r)s 


Elf  I  y,- 1  Q 

Pr\  2  a,  2  a,) 

(51) 


which  is  solved  numerically  with  the  Newton-Raphson 
method.  The  compatibility  conditions  are  once  again  found 
by  noting  that  across  the  expansion  waves  0  ^  Pf/p,  ^  1  and 
0^  Pfip,  ^  1,  which,  using  Eqs.  (50),  give 

0<u,-«,<2f-^  +  -^)  (52) 

Vy,-1  V/-1/ 

and 


<r^0. 


aen.w  the  esptasvusrm  wave 

tv  ^-txnpwrihshis  cowdvriosM  art  found  by  notnig  that 
icrvvas  the  espwntswHs  wave  L  p.  p  am)  acrowi  the 
shock  P.  p  p  >  ^  I  smg  Fqs  and  f46l  the  Mkcrwvng 
isflarwvns  are  Si’und.  aher  some  akpefwa. 


a 

•  * 

ansf 


0  s;  » 

Hw'  I  -s4t«c4.  ^-rriMwi'asm  nee  TVs  «»*  rs  eeactfy  the 
same  as  cam  t  n »  wtih  the  rrawihufmarava  f  £ 


3  2  Acotalic  Approximation 

The  aolution  to  the  Riemann  proMem  becomes  easier  to 
ohum  when  the  initial  conditions  are  such  that  the  waves  R 
and  L  shown  in  Fig.  7  are  to  weak  that  linear  acoustic 
theory  can  be  used  This  happens  when  the  distance,  in 
scene  sense,  between  the  two  constant  states  r  and  /  is  small. 
The  waves  R  and  L  can  then  be  treated  as  acoustic  waves 
with  the  pressure  ratios  acron  them  given  by  the  simple 
retatsosK. 


P<-  P-*  ^  7rp,p.ifr-it,).  (54a) 

Pf- P- (54b) 

r  cesdsNMSf  Eqs  <  54a  t  and  ( .54b  L  one  finds 


=  -  p.^  ^  ~,.p,p,a. 

'  P'P'»A'iy,  '.PrP.*~Jy,p,p^  (55) 
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The  densities  behind  these  waves  are  given  by 


where 


P?  =  (V,-l)/(r,  +  l),  i  =  rJ.  (57) 


3.3.  Reacting  Mixture  of  Caloricaliy  Perfect  Gases 


So  far,  the  classical  Riemann  problem  for  non-reacting 
inviscid  flow  has  been  considered.  The  solution  to  this 
problem,  as  mentioned  before,  is  a  self-similar  solution.  i.c.. 
depending  on  xjt  only.  For  the  case  of  a  simple  reacting 
mixture  the  nonlinear  system  of  equations,  that  needs  to  be 
solved,  is  of  the  form 


dU 

dt 


+ 


8F{U) 

dx 


=  G(f/), 


(58) 


where 


(  ^  \ 

{  f  \ 

u= 

pu 

p(c+;«^)  1 

,  F(C/)  = 

P“  +  P  1 

pu(e  +  ^  u^)  +  pu  1 

^  P.’  ) 

[  pzu  ) 

and 


(59) 


G(U)  = 


(60) 


This  is  written  using  the  Eulerian  formulation,  but  one 
obtains  a  system  of  exactly  the  same  form,  if  the  Lagrangian 
formulation  is  used.  The  Riemann  problem  solution. 


FIG.  8.  This  is  a  typical  wave  pattern  resulting  from  the  breakup  of 
the  initial  discontinuity  of  the  Riemann  problem  for  the  case  of  a  simple 
reacting  mixture.  The  solution  is  no  longer  self-similar. 


dewnbed  i*i.  «  tiM  it»e  tw.*  i‘r«t't4ttj  cMm  i  t  c.i  t  >  -  ii 
For  the  reat.'iing  tat*  t*s  the  til  umIhuui 

gas  are  allotaexJ  only  aiit**.*  mstiisii’i 
not  across  shocks  The  kt4utKS»i  ii.'  thi»  fttutstem.  is  «ii«e 
complicatcxJ  atsd  no  ksAgrt  sirfl  In  ♦  tf  F  »  tvpnal 

wave  pattern  ts  shown  The  sht«k  and  rs|VMJfcit»ti  waves  a>e 
curved  in  the  (i,  tt  piane  »e  they  ate  ao«tle!tat»tt*  th* 
solution  to  this  fmetaiusd  Riewann  has  Howi 

worked  out  by  Matanu  Ben  Anri  1  !  ]  It  »  show-ti  that  tin 
solution  approaches  the  soluison  tii  thr  ciaaattwi  Rtomartti 
problem  for  the  non-rracun*  caae  m  the  htnit  •  ci  and 
t-*0.  The  more  cxsmpiscaied  gnseralmod  Rsnmartti  aolvn 
given  by  Ben-Art/i  provides  higher  <*fdrf  at*  utat  y  ovn  the 
usual  non-rcacting  solver  Numerical  cspr»»*»e«ii  were  pet 
formed  using  the  preseni  adaptive  l.aftat^pat>  athemt  w'lth 
both  Riemann  solvers  It  was  lound  iHai  ihe  simplrt  stilvrt 
gave  results  which  were  just  as  good  The  at**^ai«»ei  td  ilw 
various  waves  was  captured  nunserx^Uy  quite  affutateli 


4.  NI  MKRICAI  RtNI  IN 

A  number  of  test  cases  were  riin  using  this  numerrscal 
scheme.  The  cases  were  chosen  pnmanly  to  salidair  the 
code  and  to  demonstrate  its  potential  (or  solving  1{) 
problems  with  complicated  wave  interactions  The  scheme 
is  particularly  useful  for  computing  unsteady  reading  (lows 
involving  detonation  waves  and  their  inieradions 

4.1.  Sod's  Shock-Tube  Problem 

The  first  case  is  the  classical  shock-lubc  problem  It 
is  an  important  validation  run  for  the  code  The  initial 
conditions  used  are  those  proposed  by  Sod  [13]  At  iiiik 
r  =  0  a  diaphragm  at  the  location  x  =  0.5  separates  the  two 
constant  states 

p,=  1.0,  p,=  1.0,  M,  =  0.0,  .t<0.5 

(61 ) 

p,  =  0.1,  =  0.125,  «,  =  0.0.  .v>0.5. 

for  a  perfect  gas  with  y  =  1 .40.  ^=150  computational  cells 
are  used  in  this  calculation.  In  ail  the  results  presented, 
the  solutions  are  given  as  functions  of  the  Eulerian 
space  variable  x,  even  though  the  calculation  is  done  in 
Lagrangian  space.  The  Lagrangian  aspect  of  the  scheme  is 
evident  by  the  increased  density  of  computational  points  in 
compression  regions.  The  comparison  between  the  numeri¬ 
cal  solution  and  the  exact  solution  shown  in  Fig.  9  is 
excellent.  Note  that  the  expansion  wave  is  computed  with 
the  accuracy  of  typical  shock-capturing  schemes,  since  no 
effort  is  made  to  track  expansion  waves.  The  shock  wave 
and  the  contact  discontinuity  are  computed  with  no 
smearing. 
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FIG.  F.  (a>  VflioDfv  proMr  at  mmc  r»0JD  tor  Sod's  sbodi  tobc 
proNtem  <ntS  ,¥  •  1 10  Miwi(pom»«mal  oRs  (h)  Frraift  proCIt  al  tHNC 
I  ^  0  JO  ft>r  Sod'i  shock  mSt.  proMtNi  «rNll  W  «  I  $0  coMOoUtfOMl  oefls 
(c)  Denote  prodke  af  him  r>0J0  far  Sod's  shock  lake  psoMtin  snth 
¥•130  compotatioMl  ccfli  Thr  whd  hoes  ftprcjMW  the  exact  sofabom 
and  the  boxes  tepreaent  the  wawerscal  sofatsow 


In  order  to  demonstrate  the  ability  of  the  scheme  to  com¬ 
pute  complicated  wave  interactions  accurately,  the  shock 
tube  problem  is  carried  a  step  further.  Reflecting  walls  are 
assumed  present  at  the  locations  .r  =  0  and  .v  =  I.  The  com¬ 
putation  is  continued  to  see  how  the  multiple  reflections  of 
the  shock  from  the  walls  and  its  collisions  with  the  contact 
discontinuity  are  calculated.  In  Fig.  10,  the  solution  is 
shown  after  the  first  reflection  of  the  shock  wave  from  the 
wall  at  .v=l.  which  occurs  at  f  =  0.285.  In  Fig.  11,  the 
reflected  shock  has  collided  with  the  contact  discontinuity 
and  a  new  shock  wave  has  been  generated.  The  solution  at 
a  later  time  is  shown  in  Figs.  12.  The  computation  was 
carried  out  until  time  t  =  7.88.  That  corresponds  to  many 
reflections  of  the  original  shock.  In  a  real  exp)eriment  viscous 
effects  would  have  made  the  process  die  down  much  sooner. 
In  Fig.  13,  the  entropy  of  the  system  is  shown  as  a  function 
of  time.  The  entropy  is  defined  by 

s  =  \n{p/p'‘).  (62) 


?  4  6  8  1  0 


‘  b 


C  ?  4  6  6  1  0 


FIG.  IS.  (■>  PicHUR  profile  at  time  f-0.40.  (b)  Deiuiiy  profile  at 
iHBe  (•0  40.  The  ihocfc  haa  reflected  from  the  wall  at  jr  •  I . 
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0  .2  <1  6  8  10 


FiG.  It.  (a)  Pressure  profile  at  time  r  =  0.4S.  (b)  Density  profile  at 
time  r  =  0.45.  The  reflected  shock  has  collided  with  the  contact  discon¬ 
tinuity.  A  secondary  shock  has  been  generated. 


As  t-*  00,  the  system  approaches  the  state  predicted  by 
thermodynamics,  since  the  scheme  is  fully  conservative.  Any 
scheme  which  conserves  total  mass  and  energy  will  give  the 
correct  Hnal  entropy.  In  this  case  it  is  j  =  0.1168  in  the 
appropriate  dimensionless  units.  This  is  an  important  point 
worth  repeating  here.  The  conservative  character  of  the 
scheme  is  not  compromised  by  the  use  of  the  adaptive  grid 
technique. 


4.2.  Strong  Shock  IVave  Problem 

The  strong  shock  wave  problem  used  by  Woodward  and 
Colella  [18]  is  computed  with  the  present  scheme.  This 
problem  is  a  good  test  case  because  of  the  strong  interacting 
discontinuous  waves.  The  initial  condition  is  that  of  a  gas 


0  ! . . 

C  <■  »  *  • 


FIG.  12.  (•)  Prewure  proAk  M  tmm  ’•Oh]  flit  Dtwwn  iwaWt  *1 
tii»e  f  •  041 .  An  teovauc  wavt  hM  bm  ptwwwJ  tram  dw  ctdlMto*  nf  dw 
lecondary  shock  with  the  oodtact  dwcoiwwwwr;' 


FIG.  13.  Entropy  i  ••  Intp/p')  of  the  sysiein  ns  •  hnciian  of  me  for 
the  shock  tube  proMem  with  multiple  fiRectiods. 
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*  *.  »  10 


flG  Ml  («)  VtfcaWty  >»olii  wwl  (»)  ilwiir)  prolte  ii  tiai*  >-0.030 
Ibr  dk*  «rnitf  tiindl  yruftlwi 


with  specific  heal  laiio  y»  14  at  real  in  the  tube0<x<  I. 
The  niittal  demity  is  p  «  I  and  the  ptesaure  is 

^=1000.  jr<0J. 

^-001.  0.l<x<0.9.  (63) 

100.  0,9<jt<  1. 

The  resufta  are  shown  in  Fip.  14  and  IS  for  the  times 
r»  0.030  and  r»  0.038.  re^wctively:  800  computational 
ce)h  were  used  for  this  calculation. 


« 


FIG.  IS.  (•)  Velocity  profile  and  (b)  density  prifile  at  time  i  =  0.038  for 
the  strong  shock  wave  problem. 


4.3.  ZND  Detonation  Waves 

Another  test  case  is  that  of  a  steady  detonation  wave. 
The  well-known  ZND  theory  (Zel’dovich-VonNeumann- 
Doering)  for  a  steady  detonation  is  used  to  compare  with 
the  numerical  solution  obtained  using  this  scheme.  As  a  first 
check,  the  profile  of  a  steady  detonation  wave,  computed 
using  the  ZND  theory,  is  given  as  the  initial  condition  to  the 
unsteady  code.  The  solution  after  time  t  =  10  (10,000  time 
steps)  is  then  superimposed  on  the  ZND  solution  and 
compared.  The  comparison,  shown  in  Fig.  16,  is  excellent. 


FIG.  16.  (a)  Velocity  profile  at  10.  (b)  Pressure  profile  at  10.  (c)  Density  profile  at  r*  10  (6)  Towferwrr  prtHSif  f »  lb  tlir  «nM1  Imi^ 
are  the  solutions  given  by  the  ZND  theory.  The  boxes  are  the  numehcat  solutions. 


The  standard  Arrhenius  law,  given  by  Eq.  (la),  is  used  for 
the  chemical  reaction  rate  with  a  =  0,  i.e., 

i=  -KzT’np{-E/R^T). 

The  parameters  used  for  this  test  run  are 

y=1.2,  qJR,To  =  SO,  EIR,To  =  AO, 

where  the  subscript  zero  denotes  the  undisturbed  region 
into  which  the  detonation  propagates.  This  steady  detona¬ 
tion  wave  corresponds  to  an  overdrive  factor  of/>  1.6.  The 
overdrive  factor  is  defined  by 

where  D  is  the  detonation  wave  speed  and  Drj  is  the  detona¬ 


tion  speed  corresponding  to  the  Chapman  Joiifttei  pmnt 
For  dmils  on  the  ZND  theory  aee  the  book  by  Ficketi  sod 
Davis  [4]. 

The  case  of  unsteady  detonation  «ra«ct  «stl  now  be  oon 
skfered.  For  the  foDowing  cases  the  swwpliiied  Arrhenno 
chemical  rate  law  is  <«ted  (Fq  ( lb)l  te, 

-KiHiT-  T,l 

where  T,  is  a  critical  temperature  above  whsoh  the  dhenncai 
reaaioo  begins.  Figure  1 7  shows  the  evolatton  of  an 
unsteady  detorution  propagating  ta  an  Msdistarbed  region 
It  is  the  well-known  piston  problem  The  motion  of  the 
piston,  starting  at  x^O.  generates  a  rfiock  which  raiaes 
the  temperature  of  the  gas  above  the  critical  valwe 
The  chemical  reaction  begins  and  the  detonation  srave 
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cpndttmm  «•  •  wwwth  preienwe  dmfn9mtwm  «t(h  aero  ll•f(l•i 
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very  well  by  this  scheme.  N  —  200  computational  cells  were 
used.  The  numerical  values  used  in  this  run  are 

7=1.4,  9o  =  20,  re=1.2,  K=l, 

where  all  quantities  are  normalized  appropriately.  The 
normalization  is  such  that  the  initial  temperature  of  the  gas 
at  rest  is  To  =  1  or,  equivalently,  po  =  Po  in  dimensionless 
units.  The  initial  pressure  distribution  is  given  by 

Po(x )  =  0. 1 0  +  3.0  exp(  -  2Sx- ). 


5.  CONCLUSIONS 

An  adaptive  numerical  scheme  has  been  presented  for 
the  computation  of  flows  with  complicated  interactions  of 
discontinuous  waves.  Its  accuracy  and  robustness,  as 
demonstrated  by  numerical  experiments  make  it  a  valuable 
tool,  especially  for  the  study  of  unsteady  reacting  flows  with 
detonation  waves.  The  conservative  formulation  gives  the 
method  all  the  advantages  of  higher-order  shock-capturing 
schemes  and  its  adaptive  characteristic  allows  for  good 
accuracy  near  shocks  with  no  smearing  effect.  The  advan¬ 
tages  of  the  conservative  shock -capturing  schemes  are  com¬ 
bined  with  the  advantages  of  the  front-tracking  methods 
very  well  to  give  a  useful  computational  scheme. 

The  drawback  is  that  the  extension  of  this  scheme  to 
multidimensional  flows  is  not  straightforward.  The  main 
idea  of  the  scheme  is  the  conservative  front-tracking  of 
shocks  and  contact  discontinuities  on  a  Lagrangian  grid. 
The  Lagrangian  aspect  of  the  method  is  the  most  difficult 


to  extend.  The  conservative  front-tracking  aspect  can  be 
extended  and  work  in  this  area  is  in  progress. 
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Abstract 

This  paper  focuses  on  the  correlation  of  two  successive  scalar  images  for  the 
purpose  of  measuring  imaged  fluid  motions.  A  method  is  presented  for  deforming, 
or  transforming,  one  image  to  another.  Taylor  series  expansions  of  the  Lagrangian 
displacement  field  are  used,  in  conjunction  with  an  integral  form  of  the  equations 
of  motion,  to  approximate  this  transformation.  The  proposed  method  locally  corre¬ 
lates  images  for  displacements,  rotations,  deformations,  and  higher  order  displace¬ 
ment  gradient  fields,  and  applies  a  global  minimization  procedure  to  insure  a  global 
consistency  in  the  results.  An  integral  form  of  the  equations  of  motion  is  employed 
and,  as  a  consequence,  no  spatial  or  temporal  differentiation  of  the  image  data  is 
required  in  estimating  the  displacement  field.  Successive  two-dimensional  digital 
CCD  images  of  fiuid  motion  marked  with  dye,  are  used  to  verify  the  capabilities  of 
the  method.  The  utility  of  the  method  is  also  illustrated  using  a  pair  of  Voyager  2 
images  of  Jupiter. 
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1.  Introduction 


The  application  of  photographic,  CCD,  and  other  forms  of  imaging  for  the 
purpose  of  estimating  flow  velocities,  has  been  investigated  by  many  researchers  in 
fields  ranging  from  fluid  mechanics  to  vision  research.  In  the  most  common  methods 
for  measuring  fluid  flow  velocities,  the  flow  is  seeded  with  particles,  or  markers,  that 
can  be  easily  imaged  and  tracked.  An  extensive  review  of  methods  using  particle 
and  speckle  images  for  fluid  flow  measurement  is  presented  by  Adrian  (1991).  The 
estimation  of  the  motion  and  deformation  of  solids  is  closely  related  to  that  of  fluids. 
A  method  of  determining  displacements  and  stress  intensity  factors  in  solids,  using 
white  light  speckle  images  and  image  correlation  techniques,  is  presented  in  McNeill 
et  al.  1987.  In  the  absence  of  particles,  flows  have  also  been  tagged  with  a  line 
or  grid,  e.g.,  using  laser-induced  photochemical  reactions  (Falco  &  Chu  1987),  or 
laser-induced  fluorescence  (Miles  et  al.  1989).  When  this  is  not  possible,  one  can 
use  markers  that  occtu*  naturally  in  the  flow,  e.g.,  Bindschadler  and  Scambos  ( 1991) 
have  correlated  the  translation  of  distinct  surface  features  in  ice  flows  to  determine 
flow  velocities. 

Determining  motion  from  successive  images  is  also  of  interest  in  animation,  as 
well  as  the  study  of  biological  and  robotic  vision.  Most  investigations  along  these 
lines  have  taken  the  form  of  extracting  the  motion  of  objects  in  an  image  and,  as  a 
consequence,  they  focus  on  the  motion  of  rigid  objects  and  their  representations.  See 
Hildreth  &  Koch  (1987)  for  a  review  and  Miuray  &  Buxton  (1990).  This  approach 
is  somewhat  difierent  from  the  interests  of  Fluid  Mechanics  where  the  object  of 
interest  is  a  fluid,  highly  deformable  and  often  compressible.  Nevertheless,  many 
results  from  object  motion  research  apply  directly  to  the  motion  of  fluids  emd  solids. 

The  proposed  Image  Correlation  Velocimetry  (ICV)  method  that  will  be  de¬ 
veloped  in  the  present  discussion  has  roots  in  both  the  correlation  methods  used 
in  measuring  fluid  flow  and  the  deformation  of  solids,  outlined  in  Sec.  1.1,  and  the 
gradient  methods  used  in  measuring  optical  flows,  outlined  in  Sec.  1.2. 


1.1  Correlation  methods 

Several  techniques  for  determining  fluid  flow  velocities  from  particle  image  pairs 
(e.g.,  Willert  &  Gharib  1991)  employ  an  optimization  that  relies  on  some  form  of  a 
cross-correlation  function,  e.g., 


max  /  Eo(x)  Fi(()  d*x  , 
•  Ja 


(1) 


with 


C  =  x  +  a  , 


(2) 
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where  a  is  a  vector  parameter  to  be  determined  by  the  optimization  procedure  and 
A  is  the  correlation  region.  The  distribution  of  the  image  irradiance,  E(x,  t),  is 
known  at  times  to  ti,  t.e., 

Eo{x)  =  E{x,to)  and  i;i(x)  =  f?(x,ti)  .  (3) 


The  average  velocity,  ,  within  the  correlation  region  is  then  approximated  by, 


u 


A 


a 

tl  —  to 


(4) 


The  drawback  of  this  method,  having  only  two  parameters  to  quantify  the  motion, 
is  that  it  cannot  resolve  displacements  properly  where  there  are  large  displacement 
gradients  within  the  correlation  region.  Anticipating  this  problem,  and  being  very 
interested  in  displacement  gradients,  researchers  in  Solid  Mechanics,  apply  tech¬ 
niques  which  include  higher  order  deformations  of  the  displacement  field  within  a 
correlation  volume.  For  example,  McNeill  et  al.  (1987)  describe  a  method  whereby 
a  model  of  the  image  displacement  field  (mapping)  is  used  in  a  least  squares  opti¬ 
mization  procedure,  t.e.. 


min  f  [Eq{x)  -  Ei{^)]^  d^x  .  (5) 

The  affine  mapping, 

^  =  X -H  a -b  (Va)  •  dx  ,  (6) 

is  used  as  an  example  of  such  a  function,  and  the  displacement  a  and  the  four 
components  of  Va  are  treated  as  parameters  to  be  determined  by  the  optimization 
procedure.  However,  any  physically  motivated  mapping  can  be  xised  in  place  of 
Eq.  6. 


In  both  these  methods,  the  image  data  are  integrated  over  a  region  and  require 
no  spatial  differentiation.  Since,  for  two-dimensional  images,  only  a  few  parameters 
are  extracted  from  the  optimization,  these  methods  are  relatively  immune  to  noise 
and  lend  themselves  to  fast  solutions. 


1.2  Gradient  methods 


A  method  for  determining  the  velocities  of  visual  features  in  an  image  was 
presented  by  Horn  &  Schimck  (1981).  This  visual  velocity  is  termed  “optical  flow” 
to  differentiate  it  from  the  velocities  of  (material)  objects  in  an  image,  e.g.,  a  shadow 
moving  across  the  ground  has  a  perceived  velocity  that  is  markedly  different  from 
that  of  the  groimd,  and  a  rotating  featureless  disk  will  have  no  visual  velocity  at 
all.  The  fundamental  equation  iised  by  Horn  &  Schunch  to  determine  the  optical 
flow  was, 


-f-  u  •  Vf?  =  0 
at 


(7) 
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where  E  is  the  image  irradiance  (a  scalar)  and  u  is  the  optical  flow  velocity.  The 
differential  terms,  dEfdt  and  V£,  can  be  estimated  from  the  image  data  and  the 
component  of  u  along  is  calculated  using  Eq.  7.  Methods  employing  equations 
of  this  type  are  called  gradient  schemes.  Note  that  no  velocity  can  be  calculated 
using  Eq.  7  if  there  are  no  featmes  or  gradients  in  the  image,  t.e.,  if  VE  =  0. 
In  addition,  because  Eq.  7  employs  only  the  component  of  u  along  VE,  velocity 
components  along  the  equi-scalar  contours  of  E  cannot  be  determined  using  Eq.  7 
alone.  This  limitation  was  designated  the  “aperture  problem”  by  Wallach  (1976). 
The  terminology  is  somewhat  misleading  and  is  only  used  here  in  reference  to  the 
convention.  More  appropriate  names  might  be  the  “characteristic  problem,”  or  the 
“direction  problem,”  because  the  problem  is  finding  the  velocity  along  the  charac¬ 
teristic  direction,  E  =  constant.  Gradient  schemes  also  have  the  problem  that  finite 
difference  approximations  of  the  spatial  and  temporal  derivatives  are  necessary.  A 
problem  with  such  approximations  for  the  derivatives  is  related  to  the  Nyquist  sam¬ 
pling  criterion,  where  aliasing  in  the  image  data  can  effect  the  velocity  estimates. 
To  minimize  this  problem  in  taking  the  gradient,  the  motion  between  images  should 
be  less  than  half  the  smallest  local  spatial  scale,  Xe,  of  the  E— field,  t.e.. 


|u|  ih-tp)  1 

Xe  2 


(8) 


(c/.  Eq.  7),  where  {ti  —  to)  is  the  time  between  images. 


The  uncertainty  of  the  so-called  aperture  problem  can  be  solved  in  some  cases 
by  applying  constraints  to  the  motion,  e.g.,  the  motion  is  of  a  rigid  body  (see 
Murray  &  Buxton  1990,  for  example),  or  a  limited  class  of  deformable  bodies  (see 
Terzopoulos  &  Metaxas  1991). 


Horn  &  Schunck  (1981)  applied  a  global  constraint  to  Eq.  7  (in  two  dimensions), 
by  solving  for  u(x,t)  using  an  optimization,  *.e.. 


min 

u(x,t) 


(9) 


where  a  represents  the  constraint  cost  function  in  the  optimization  process,  and 
k  balances  the  relative  cost  of  <t  and  Ekj.  7.  In  particular,  Horn  &  Schimck  chose 
smoothness  as  a  constraint,  t.e.. 


(10) 


The  idea  of  including  constraints  in  the  optimization  process  that  determines  the 
velocity  field,  over  an  area,  is  important  in  the  context  of  the  method  to  be  discussed 
below.  Note  that  the  constraint  in  Eq.  9  need  not  be  included  in  the  optimization 
integral.  Instead,  it  could  be  included  as  a  feature  of  the  optimization  technique. 
See  Murray  &  Buxton  (1990). 
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2.  PropoMd  metbodologjr 

A  fucceMMMi  of  tnumai  can  rtprmmi  anytiiiof  from  the  moiAoo  of  cam  oo  a 
highway  to  the  tranaiMirt  of  a  dye  outfhcf  to  wausr.  We  take  ibe  view  that  givoo 
succeaaive  image  rcprcaentationa,  there  exiata  a  iranaCormaUoo.  or  mappiAg.  of  the 
local  image  intensity  data  that  takes  ooe  image  to  the  nest  In  many  rasas,  while 
the  equ^ion  of  motMMi  of  the  imaged  field  nuty  be  known,  the  ma{>timg  taking  one 
image  to  another  may  not  be.  Successive  images  confined  with  the  oquattons  of 
motion,  however,  often  allow  us  to  approximate  the  mapping 

The  mapping  of  one  image  to  the  next  can  be  developed  by  coos^ermg  the 
Lagrangian  displacement  field  ((x,()  of  the  image  sequence.  Specifically,  if  x  is 
the  coordinate  of  a  point  on  an  image  at  some  initial  time  ((x.f)  represenis 
the  coordinate  of  this  point  in  a  subsequent  image  recorded  at  a  later  time  ( .  If  we 
imagine  the  image  sequence  as  the  result  of  a  continuous  recording  process,  we  can 
assign  a  Lagrangian  tma^e-flow  velocity  field,  referred  to  as  “optical  flow*  in  the 
discussion  and  literature  cited  above,  t,e., 


u[^(x,t),t]  =r  ,  (11) 

to  the  continuous  displacement  field  ((x,t)  that  takes  an  initial  point  x  in  the 
image  recorded  at  time  to ,  to  the  point  ^(x,  t)  on  the  image  recorded  at  time  t . 
We  recognize  that,  for  the  case  where  the  images  represent  fluid  flow,  e.g.,  successive 
images  of  a  convected  scalar,  the  image-flow  velocity  field,  u(x,  t) ,  may  be  quite 
difierent  from  the  fluidrRow  velocity  field  Uf(x,t).  The  extent  to  which  the  former 
represents  a  good  approximation  for  the  latter  is  a  separate  issue  that  can  only  be 
addressed  in  the  context  of  the  details  of  the  particular  imaging  process  and  the 
fluid-flow  field. 

In  the  proposed  implementation,  local  series  approximations  for  the  displace¬ 
ment  mapping  are  used  in  conjunction  with  an  integral  form  of  the  equations  of 
motion.  A  global  nonlinear  correlation  (optimization)  process  is  employed  to  esti¬ 
mate  the  image-flow  velocity,  vorticity,  deformation  rate,  etc.,  of  the  imaged  data 
field.  “Series,”  in  this  discussion,  will  denote  “Taylor  series.” 

In  the  context  of  fluid  mechanics  measurements,  we  will  focus  on  images  of 
continuous,  passive,  convected  scalars,  e.g.,  dye  markers,  carried  by  a  fluid.  As  will 
be  illustrated  using  a  pair  of  Voyager  2  images  of  Jupiter  (Sec.  5),  however,  any 
marker  in  the  flow  can  be  used. 

The  method  will  be  developed  for  three  dimensions  and  can  yield  three- 
dimensional  velocity  fields.  The  method  can  also  obtain  two-dimensional  velocity 
fields  from  images  of  two-dimensional  flows.  In  a  concession  to  the  limitations  of 
typical  data  acquisition  systems  today,  however,  the  method  will  be  applied  here 
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to  a  two-dimensional  flow  and  also  to  two-dimensional  slices  of  three-dimensional 
flows.  We  note  that  in  some  cases,  two-dimensional  imaging  devices  can  be  used  to 
obtain  approximations  to  three-dimensional  image  data,  e.y.,  Dahm,  et  al.  1991.  A 
short  discussion  of  the  implications  of  correlating  two-dimensional  slices  of  three- 
dimensional  data  is  presented  in  Sec.  4. 


2.1  Fluid  displacement  and  equations  of  motion 

To  see  bow  the  image-flow  velocity  field  can  be  calculated  from  three- 
dimensional  image  data  sets,  spaced  in  time,  first  consider  a  Lagrangian  description 
of  a  flow  being  imaged.  Figure  1  illustrates  the  motion  of  fluid  particles  within  a 
volume,  V .  Fluid  elements  at  x ,  in  a  neighborhood  V ,  at  time  to  <  are  convected  to 
locations  ^(x.  ()  at  a  later  time  t.  The  displacement  field,  ((x,t),  can  be  thought 
of  as  a  transfmination  of  the  field  x,  id  time  to,  to  ((x,t).  Given  the  image 
displacement  field  C(x,t).  the  image-flow  velocity  field  is  then  given  by  Eq.  11. 


Fig.  1  Motion  of  a  fluid  volume. 


Using  this  Lagrangian  field,  ((x,  t) ,  one  could,  in  principle,  integrate  the  equa¬ 
tion  of  motion  of  the  imaged  scalar,  t.e., 

^  u  •  Vc  =  PV*c  ,  (12) 

at 


ci[^(x,ti)]  -q)[^(x,to)]  -V 


V*c[^(x,t),t]  dr  =  0  , 


(13) 


to  obtain 
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where 

Oo(x)  s  c(x,to)  and  ci(x)  s  c(x,ti)  (14) 

represent  the  c(x,  t)‘field  at  times  to  and  ti ,  respectively,  and  V  is  the  expropriate 
diffusion  coeflBdent. 


In  the  first  two  examples,  the  motion  of  food  coloring  in  glycerine  (Sec.  3)  and 
dilute  fiuorescein  in  water  (Sec.  4),  are  examined.  In  th^  fiows,  the  diffusion  of 
the  dye  meurkers,  in  the  time  interval  between  successive  images,  is  relatively  small 
and  may  be  neglected,  t.e.. 


(<i-to)P 


<  1  , 


(15) 


where,  ti  —  to  =  0.1  s  is  the  time  between  images,  /  =  50  to  80  pm  is  the  imaging 
resolution,  and  the  diffusion  coefficients  are  no  larger  than  V  —  10~*m^/8.  In 
addition,  we  note  that  the  Schmidt  number  is  large  in  both  fiows,  t.e.. 


Sc  S  iz/B  >  lO’  , 


(16) 


where  u  is  the  kinematic  viscosity.  In  the  first  example  of  the  dye  marker  in  glyc¬ 
erine,  the  fiuid  fiow  is  two-dimensional,  as  is  the  image,  and  the  image-fiow  velocity 
field,  u(x,t),  is  a  good  representation  of  the  fluid-flow  velocity  field,  Uf(x,t).  In 
the  second  example,  both  the  flow  velocity  and  the  imaged  scalar  field  are  three- 
dimensional,  while  the  image  is  two-dimensional.  As  we  will  discuss,  the  image-flow 
field  need  not  necessarily  represent  the  fiow  velocity  field,  in  that  case.  In  the  third 
example,  the  motion  of  the  imaged  quantity  in  the  Jovian  atmosphere  (Sec.  5)  does 
not  follow  any  simple  equation  of  motion.  In  that  example,  the  derived  image-flow 
velocity  field  can  be  expected  to  be  an  even  poorer  representation  of  the  fluid-flow 
velocity  field. 


In  cases  where  the  diffusion  of  the  imaged  scalar  can  be  ignored,  Eq.  13  becomes 

ci[€(x,ti)] -co[((x,to)]  =  0  .  (17) 

Equation  17  represents  a  significant  simplification  over  Eq.  12,  its  differential  coun¬ 
terpart.  It  contains  no  spatial,  or  temporal,  derivatives  and  suffers  few  of  the 
drawbacks  associated  with  the  gradient  methods  discussed  earlier  (Sec.  1.2). 


Using  the  integral  equation  of  motion  (Eq.  17)  in  place  of  the  differential  equa¬ 
tion  of  motion  (Eq.  7)  in  the  optimization  (Eq.  9),  and  generalizing  the  optimization 
to  three  dimensions  then  yields  an  expression  for  determining  ((x,to)  ((x,ti), 

t.e., 


min  f  (  {  ci[C(x,t|)] --co[^(x,to)]  }*  +  d®x.  (18) 
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In  the  spirit  of  the  correlation  methods  discussed  in  Sec.  1.1,  where  the  type 
of  motion  within  the  correlation  volume  is  limited  to  translation  alone,  the  present 
method  restricts  ((x,  to)  and  ((x,  ti)  to  the  first  few  terms  of  a  series  approximation 
for  ((x,  t) .  While  many  representations  of  the  displacement  field  could  be  employed, 
a  (Taylor)  series  representation  is  used  because  the  first  two  orders  in  the  series 
expansion  correspond  to  physical  fluid  mechanical  quantities,  t.e.,  the  velocity  vector 
and  the  velocity  gradient  tensor.  More  importantly,  the  series  approximation  has 
the  additional  benefit  of  enforcing  smoothness  in  the  displacement  and  displacement 
gradient  fields  within  a  correlation  volume. 


2.2  Displacement  field  and  kinematic  quantities 


In  the  case  of  fluid-flow  images,  the  quantity  ((x,  t)  is  a  complicated  nonlinear 
function  of  the  imaging  process,  the  nonlinear  convection  dynamics,  and  x.  Local 
estimates  of  this  function  can  be  made  by  Taylor  series,  expanding  ((x,  t)  in  space, 
at  some  time  t,  in  an  image  correlation  volume,  V.  This  yields, 

^(x, t)  -  ^(xc-, t)  +  {x-Xc)-  V^(xc; t)  +  ^  ((x  - Xc)  •  Vf  +  (19) 

In  this  expression,  x<;  denotes  the  center  of  the  image  correlation  volume,  V,  at 
time  <0,  and  V((Xc;t)  denotes  the  gradient  of  ((x,t)  with  respect  to  x,  evaluated 
at  Xc.  Figure  2  plots  the  number  of  parameters  used  in  the  optimization  process  as 
a  function  of  the  order  used  in  the  series  expansion,  for  two  and  three  dimensions. 
Figure  3  illustrates  the  effect  of  the  various  orders  of  the  expansion  on  a  two- 
dimensional  square  “volume.” 

Using  a  finite  difference  approximation  in  time  for  the  velocity,  Eq.  11,  and  the 
series  representation,  Eq.  19,  evaluated  at  times  to  ^<1  yields  an  estimate  for 
the  velocity  within  the  correlation  voliune,  t.c.. 


(20) 


where  to  <  t  <  ti.  Similarly,  taking  the  spatial  gradient  of  Eq.  20  yields  an  expres¬ 
sion  for  the  velocity  gradient  tensor  within  the  correlation  volume,  t.e., 

Vu[«x;(),t]  s  ^u[«(x;t),(]  =  — -  ■  (21) 


Vorticity,  divergence,  and  strain  rate  can  then  be  obtained  from  the  components  of 
the  estimated  velocity  gradient  tensor,  Eq.  21. 


ber  of  parameters 
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Order  of  series  expansion 


Fig.  2  This  plot  shows  the  rt^id  increase  in  the  number  of  parameters  used  in  the 
optinuzation  procedure  with  increasing  order  of  the  series  expansion.  is 
for  a  2-D  expansion  and^o”  for  3-D. 

There  is  some  freedom  in  choosing  the  coordinate  transformation  at  the  initial 
time  tg,  {(x,  to)>  Om*  choice  is  to  have  the  coordinate  description  at  the  initial  time 
to  correspond  with  the  local  Eulerian  coordinates  at  that  time,  t.e., 

^(x,to)  =  X  .  (22) 

In  terms  of  the  series  expansion,  Eq.  19,  this  means  that 

Ve(xc;to)sI  ,  (23) 

where  I  is  the  identity  tensor,  and  all  other  higher  order  derivative  terms  are  iden¬ 
tically  zero. 


Translation  (no  deformation) 


Linear  deformation 


Quadratic  deformation 


Cubic  deformation 


Combined  translation  and  deformation 


Fig.  3  The  effect  if  translation  and  various  orders  of  deformation  on  a  two- 
dimensional  square  ‘Volume.”  See  Sec.  2.2. 
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3.3  Seeking  a  global  solution 


A  solution  for  the  coefficients  of  the  series  expSLnsion,  t.e.,  ((Xc;t),  V((Xc;t), 
etc.,  can  be  <^>tamed  in  a  neighborhood  around  Xc  using  the  expansions  for 
((x,to)  and  ((x,t)  from  Sec.  2.2.  The  unknown  coefficients  of  the  series  ex¬ 
pansion,  ((x«;t),  V((Xc;t),  etc.,  are  treated  as  parameters  in  an  optimiza¬ 
tion  process.  To  minimize  the  difference  between  two  data  sets  (images),  in 
a  least  squares  sense,  for  a  single  correlation  volume,  we  use  the  optimization, 
Eq.  18,  in  conjunction  with  the  series  approximations  developed  in  Sec.  2.2,  t.e., 

mm  f  ({ci[^(Xc;ti)-|-(x-Xc)  VC(Xc;ti)-f...]-co[x]}*-t- 


•f  ka 


d’x 


(24) 


The  optimization  implied  in  Eq.  24,  combines  many  of  the  best  features  of  corre¬ 
lation  methods  and  gradient  methods,  while  eliminating  many  of  the  deficiencies. 
Specifically,  this  optimization  method  has  high  immunity  to  noise,  uses  the  equa¬ 
tions  of  motion,  can  incorporate  constraints,  requires  no  differentiation  to  calculate 
the  displacement  field,  and  can  capture  displacement  gradients  within  a  correlation 
volume. 


In  principle,  a  single  correlation  volume  covering  the  entire  image  and  a  series 
approximation  of  a  high  enough  order  can  be  used  to  capture  the  entire  image 
displacement  field.  In  practice,  however,  employing  a  series  approximation  beyond 
the  third  order  (cubic)  term  is  impractical  because  of  the  rapid  increase  in  the 
number  of  parameters  in  the  optimization  process  with  increasing  order  (see  Fig.  2), 
and  the  associated  increase  in  the  computational  time  and  complexity.  In  the 
present  calculations,  when  the  quadratic  term  is  not  sufficient  to  capture  the  image 
deformation  over  the  entire  flow  field  using  a  single  volume,  as  is  usually  the  case, 
several  series  expansions  residing  in  smaller,  adjacent,  correlation  volumes  are  used 
in  place  of  the  single  large  volume. 


To  construct  a  global  optimization  using  a  number  of  local  series  expansions,  we 
require  that  neighboring  correlation  volumes,  with  independent  series  expansions, 
must  yield  consistent  results.  In  the  present  method,  we  use  the  expansion  for  the 
displacement  field  about  one  correlation  volume  to  estimate  those  of  its  neighbors. 
The  displacement  field  of  these  neighbors  is  also  estimated  in  terms  of  their  own  se¬ 
ries  mqpansions.  The  root-mean-square  difference  between  displacements  estimated 
by  neighboring  correlation  volumes  is  applied  as  a  constraint  cost  fimction.  Since 
it  is  necessary  to  refer  to  a  number  of  series  expansions,  it  is  useful  to  define  the 
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Fig.  4  The  kl  •  kS  solid  circles  denote  the  points  in  Vi  used  by  the  constraint  cost 
function  <Ti .  The  empty  circles  denote  their  counterparts  estimated  by  the 
neighbors  V^i  -  Vj4. 


taylor  series  for  ^(x,  t)  in  a  neighborhood  Vi  centered  about  Xc^ .  as 


^i(x;t)  =  ^(xc,;t)  +  (x-Xci)  V^(xci;t)  +  ^I(x-XcJ- V]*€(xc,;t)  +  ---  •  (25) 

When  series  expansions  about  miiltiple  points  are  then  employed,  the  minimization, 
Eq.  24,  is  modified,  <.e.. 


mm 

€(X«j;tl),V^(Xe 


/  {  {ci[C<(x;ti)]  -  co[x]  }  +  fcaM  d®x ,  V*  .  (26) 


This  minimization  is  performed  within  all  the  Vj,  simultaneously,  and  the  square 
of  the  constraint  cost  function. 


|^<(x*;0“Ci(xfc;t)p 

i  * 


(27) 


12 


is  applied  to  provide  global  continuity  of  the  solution,  denotes  the  series  expan¬ 
sion  about  the  “j”  neighbors  of  V,,  and  Xk  denotes  the  "1:”  points  of  comparison 
between  the  solutions  and  .  See  Fig.  4.  In  the  present  method,  eight  points 
about  each  correlation  volume  are  used  for  comparison,  three  with  each  of  four 
neighbors.  These  eight  points  are  sufficient  to  define  the  series  expansion  coeffi¬ 
cients  of  the  correlation  volume,  up  to  quadratic  order. 

The  present  implementation  of  the  method  can  solve  Eq.  26  for  two-dimensional 
fiow  up  to  the  cubic  term  in  the  local  series  expansions,  but  the  series  is  usually 
truncated  at  quadratic  order.  The  optimization  of  Eq.  24  is  accomplished  using  a 
multidimensional  minimization  process,  with  image  data  between  pixels  estimated 
using  bilinear  interpolation.  See  for  example  Press,  et  al.  (1988). 


2.4  Minimization  parameters  in  two  dimensions 


Typical  CCD  imaging  technologies  today  are  limited  to  two-dimensional  (spa¬ 
tial)  data.  This  is  not  a  problem  if  the  flow  being  imaged  is  also  two-dimensional. 
This  section  describes  how  the  method  is  applied  in  two  dimensions.  First,  we  de¬ 
velop  the  terms  of  the  series  expansion,  Eq.  19,  for  two-dimensional  flow.  With  the 
two-dimensional  vector 

=x-Xc  ,  (28) 

as  the  position,  x,  relative  to  the  center  of  the  correlation  volume,  Xc,  the  terms 
of  the  series  expansion  at  a  time  ti  appear  as  a  constant  term. 


^(Xc;<i)  = 


ai 

aoj 


(29) 


where  the  q<  are  the  vector  coordinates  of  the  center  of  the  correlation  volume  at 
the  time  ti ,  i.e., 

ai=^i{xc;ti)  ,  (30) 

a  linear  term, 

r  /v«  •  /v«  n  r  /Si  n 

(31) 


ttl.l  OCl,2 

Si 

02,1  02,2 . 

S2. 

where  the  Otj  represent  the  first  order  deformations  and  rotations  of  the  image 
field  within  the  correlation  volume,  t.c.. 


dXj 


(32) 


a  quadratic  term. 


0'1.22 


I _ 

• 

^1^2 

■ 

L  <^2* 

(33) 
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where 


a  cubic  term, 


otijk 


dxjext  ' 


(34) 


0=1,112 

0=1,122 

0=1.222!  . 

1 _ 

[o!2,lll 

0=2,112 

0=2,122 

0=2,222] 

L  J 

.  (35) 


with 


a  =L(  ^  \  ^^»(Xc;ti) 

-  Z\\i  +  j  +  k-3)  dxjdxkdxi  ’ 


(36) 


and  so  on  for  higher  order  terms.  The  aijt  ciijki  &re,  respectively,  related 
to  the  second  and  third  derivatives  of  the  displacement  field  within  the  correlation 
volume,  t.e.,  by  Eqs.  34  and  36. 

The  velocity  and  velocity  gradient  (Eqs.  20  and  21),  can  also  be  written  in 
terms  of  the  parameters  of  Eqs.  29  and  31  and  the  series  expansions  at  times  to 
and  tj ,  t.e., 

■■ai-Xc' 


where 


and 


-ti 


.012 -Vc 


Vu(xc,t)  = 


du/dx 

du/dy]  _  1 

0=1,1  - 1 

0=1,2 

dv/dx 

dv/dy\  ti-to 

0=2,1 

0=2,2  -  1 . 

(37) 


(38) 


(39) 


Alternatively,  the  velocity  gradient  can  be  written  in  terms  of  the  in-plane  vorticity 
and  rate-of-strain  tensor,  t.e.. 


0 

1  *xx  ^xy 

ufg/2 

0  j  ^ 

r 

Sxy  «yy . 

(40) 


where  (Vg  is  the  vorticity,  t'.e.. 


u,  = 


_  0=2.1  -  Ql.2 


tl  —  to 


(41) 


and  Sxx,  «yy,  and  Szy*  are  the  components  of  the  rate-of-strain  tensor,  t.e.. 


^  Qi.i  -  1 
tl  —  to 


(42) 
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and 


_  aig  -  1 
h-to 


1  aig  +  tt2,i 
•»"2  ti-to 


(43) 

(44) 


An  interesting  quantity  to  consider  is  the  second  invariant  of  the  rate-of-strain 
tensor  (see  Cantwell  1992,  for  example),  t.e.. 


-2q,  s  +2»J,  +  »5,  . 


(45) 


3.  Couette  flow  between  concentric  cylinders 

An  apparatus  to  generate  a  Couette  flow  between  concentric  cylinders  was 
fabricated  for  the  purpose  of  testing  the  method.  The  cylinders  were  made  from 
248  mm  lengths  of  stock  Plexiglas  tubing.  The  inner  and  outer  radii  of  the  annular 
region  between  the  cylinders  were  nominally  25.2  mm  and  40.9  mm.  The  cylinders 
were  from  stock  Plexiglas  tubing,  so  the  uncertainties  in  the  radii  were  ±lmm.  The 
outer  cylinder  was  rotated  with  a  rotation  rate  of  1.1  rad/s,  with  the  inner  cylinder 
stationary.  In  this  example,  employing  a  dye  marker  in  glycerine,  the  fluid  flow  is 
nominally  two-dimensional  and  the  marker  follows  the  flow.  Hence,  the  image-flow 
velocity  field,  u,  can  be  accepted  as  a  good  representation  of  the  fluid-flow  velocity 
field,  Uf. 


Images  were  recorded  using  a  Texas  Instruments  Multicam  MC-1134GN  Multi- 
Mode  B/W  Camera.  The  data  were  stored  digitally  using  an  in-house  multiple  frame 
grabber  (12-bit  A/D),  designed  by  Dan  Lang  and  Paul  Dimotakis  of  GALCFT,  set  to 
record  up  to  28  of  the  1134  x  468  pixel  gray  level  images  from  the  camera,  spaced 
by  100  msec  (adjtistable  between  33  and  267  msec).  Because  the  horizonted  and 
vertical  spacing  of  the  pixels  were  not  equal  on  this  CCD,  grid  spacings  and  image 
correlation  volumes  with  a  ratio  of  1:1.74  (vertical:horizontal  pixels)  were  used  to 
yield  a  uniform  spacing  of  the  data  in  the  real  image  plane.  Flow  visualization 
was  performed  by  randomly  distributing  red  food  coloring  (dye)  on  the  surface  of 
the  fluid.  To  provide  backlighting  for  the  dye  marker,  the  fluid  beneath  the  surface 
contained  a  translucent  white  suspension  of  3  fim  aluminum  oxide  ( AI2O3 )  particles 
in  glycerine.  When  illuminated  from  the  side,  this  provided  nearly  uniform  white 
backlighting  for  the  dye  being  imaged  on  the  surfsce.  Because  of  the  depth  of  field 
of  the  imaging  and  the  high  density  and  uniform  distribution  of  the  aluminum  oxide, 
scattering  from  individual  particles  was  not  detectable  in  the  video  images. 


15 


Fig.  5  Iiiitiiil  phicciiUMit.  of  s(!ri(\s  <!xp;uisiou  iK'iglihorlioods.  E;u  ii  (]riiot<‘.s  a 

control  volunK*.  Tlu;  .siiuill  circh'  at  the  (-(niter  of  (-aeli  control  volume  denotes 
the  center,  or  control  point,  of  a  seri(‘s  expansion. 


In  the  pr(!sent  investigations,  only  the  outer  cylindi-r  w;ts  rotated,  hence  the 
velocity  field  can  he  written  ;ts  (c.7..  Schlichting  1979). 


ii,{r.O)  H„{r.O)  rlr,-rjr 

— -  =  0  and  -  =  - 

r„/r,  -  r,/r,. 


(ir.) 


where  r.  f> .  a,  .  and  iin  are  the  radial  and  angular  ]>o-'itions  and  velociiii  ^  le^Hn  . 
tively.  il,,  is  the  rotation  rate  of  the  outer  cylinder,  and  r,  ;uid  r,,  aie  itn  uuiei 
and  outer  radii  of  the  cylinders  In  this  flow,  the  diveigen< c  is  zero.  1  » 

V  u  =  I)  .  (  17 1 

and  the  vortu  itv  is  iiniforni.  1  »  . 


:  1? 


^  ,{r.O] 


f 


V  •  u!  e  i 


(  l^'i 


I'k;  r>  nurnt  of  ICM)!!!*- latri  .illuwin!:  “uly  ti.m'-l.iiion 

F'>r  ilit<  ‘'iiuj'lf'  i<  *it  «  nin«'  forf<|;»iion  volum*  '*  <-<1  1>\  J'*  )tixi  1'-  vritit  ally 

,vlwl  .1  {nx*'!*-  Iiori/ontaliy.  wrrr  to  tlir  flow  'l]i>  luraliou'-  of  (li< 

imat^r  rorrt'lafion  at  tli*'  tiinr-  ;u«'  I<-<1  m  Kjj:  %  ']  hr  T*  vult‘. 

of  flir  rorrrlafif>n  jirorr<>.  allowuic  only  <l»‘«plaf rnirnt  of  ti)<  fourlatmti  volumr 
i/rroflt  orilf'f  srri<^  ox|>.ui'*iou)  shown  m  Fir  G  FtRuro  7  •frjiion'.tiatr'-  how  iIk 
flow  !•.  rapfiirorl  whrii  th*'  *  orr«’l.iJion  prorrs.s  is  rxtrn'h  rj  to  iji<  liiqlior  oiih  r 

t*'rn)s  in  ih*'  srri^s  rxp.xii'ion  To  (|ii;u)tify  tins  nnprovrjnoni .  th'’  valnr  of  th'’ 
tmtnnn/atii>n  f'lnrtional  not  inrlii'hnc  ihr  rontrihof!  o  ^  'hr  <  on'-traints.  jv  j)lottf'<l 
in  Fij;  8  ;xs  a  fniirtion  of  th*"  ofl^r  of  th*'  son'-’-  ■  .  A-  <  an  I"'  srrji,  th' 

Cr»  at*'st  improvrtj)''nf  rrali/ofl  with  th*"  mtro'l  .r  (1<  fomiatioii''  in  lli«’ 

( orr'  Iaf  ion  pro* 

Th*'  ff’siilt'  ii'^inc  ’1)*'  "Tf  lafion  iii'  tho'l  ai*  <  •  ii'-'l  with  th*’  t h'-on  t j*  al 
{ atialyti*  al  1  t wo-dini'  n'-ional  vahn''^  for  ihr  voifinfy  and  'iiv'  rc  nrr  in  Tahl**  1 
Sin*  *'  both  tlif'  vorii*  ify  and  Ihr  fhv*TC*'n*  r  arr  iimfonn  in  thr  analytiraJ  solution, 
only  a  singlr  vain*’  i>  pr*''*'nfr<l  Th*'  nnrrrtainty  in  thr  llu'on’tiral  vorfirity  is  a 
r<"'nlt  fif  th*'  rrrrnf ri(  ity  of  fh*"  rylin'h  rs  us*-*!  in  fhr  rxp'  rim'  nt  Th*'  nn<  *'rtainti*'s 
in  tlir  rx[>*"nni*'»if al  vain*"*;  arr  on*'  •itainlar*!  *l*'vialion  Ib'cansr  this  flow  has  a 
lu'arly  linrar  v*"io  ify  profil*'.  w*’  saw  only  small  rhangrs  m  thcs*'  rsfirnatrs.  h<-von<l 
th*'  hru'ar  or*l*'r  Soni''  of  th*'  “nn* *'rtainty  in  tlw  *'xp*'i inirnial  r*'snlts  for  th*' 
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I  AIU  F  1  (  onipari'on  of  i h*  or*  1 1< ,»]  .ui*!  *  xjM’riiii'  nt.il  \iir»i<)t\  .in  i  'ii\' i  >:•  n*  < 
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Order  of  series  expansion 

FlO.8  Difference  between  Couette  flow  images  under  the  mapping,  quantified  by 
the  value  of  the  optimization  functional  (arbitrary  units),  as  a  function  of 
the  number  of  terms  in  the  two-dimennonal  series  expansion. 
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4.  Cylinder  wake  flow.  Two*dimensional  slices  of  three-dimensional  data 

In  this  section,  we  present  the  resnlts  of  applying  the  two-dimensional  cor¬ 
relations  to  two-dimensional  slices  of  a  three-dimensional  flow  in  the  wake  of  an 
impulsively  started  circular  cylinder.  Here,  the  cylinder  is  1.75  cm  in  diameter  and 
45.5  cm  long.  It  is  drawn  through  a  distribution  of  fluorescein  dye  in  water  at  a 
speed  of  1.27  cm/s.  The  Reynolds  number  in  this  case  is, 

Ud 

Re  =  —  «  220  ,  (49) 

1/ 

where  U  is  the  cylinder  speed,  d  is  the  cylinder  diameter,  and  u  is  the  kinematic 
viscosity. 


Fig.  9  Initial  placement  of  series  expansion  ndghborhoods.  Each  square  denotes  a 
control  volume.  The  small  circle  at  the  center  of  each  control  volume  denotes 
the  center,  or  control  point,  of  a  series  expansion. 
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The  CCD  camera  and  data  acquisition  system  are  the  same  as  for  the  Cou- 
ette  flow  test  case  (Sec.  3).  Laser  sheet  illumination  is  provided  by  a  Continuum 
model  YG661-10  frequency-doubled  YAG  laser.  The  laser  was  operated  at  532  nm, 
300  mJ,  5  ns  pulses  width,  at  a  rate  of  10  Hz.  The  flow  here  is  three-dimensional  in 
both  the  velocity  field  and  scalar  distribution. 


Fig.  10  Displacement  of  grid  after  100  ms,  estimated  using  the  nonlinear  correlation 
process. 

Figures  9-14  demonstrate  the  method  on  images  of  a  vortical  structiure  forming 
in  the  wake  of  the  cylinder.  These  images  were  taken  after  the  cylinder  had  traveled 
about  8  diameters.  The  image  at  the  initial  time  is  shown  in  Fig.  9,  and  100  ms  later, 
in  Fig.  10.  In  this  case,  the  series  approximation  used  in  the  correlation  process  was 
expanded  to  quadratic  order.  Figure  11  shows  the  displacement  of  the  centers  of 
the  correlation  volumes. 

The  two-dimensional  vorticity  is  displayed  in  Fig.  12.  A  large  vortical  region 
can  be  seen  in  the  wake  of  the  cylinder.  The  two-dimensional  divergence,  presented 
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Fig.  11  Di.sphiccnicnt  of  centers  of  grid  over  100  ins. 


in  Fig.  13,  exposes  the  three-dimensionality  in  the  flow.  Figiire  14  plots  the  second 
invariant  of  the  rate-of-strain  ttiiisor.  Note  the  region  of  strain  (rate)  tluat  seems 
to  follow  the  periphery  of  the  large  vortical  structure.  This  could  be  a  region  of 
vorticity,  from  the  previously  shed  vortical  structure,  that  is  being  strained  around 
the  current  one. 


As  a  general  otiservation,  an  important  issue  arises  when  imaging  a  two- 
dimensional  (planar  imagi )  slice  of  a  three-dimensional  field  of  a  continuous  scalar. 
c{'K.t).  as  in  tlie  previous  example.  An  out-of-plane  component  of  the  fluid-flow  ve¬ 
locity.  Uf ,  coupled  with  an  out-of-i)lane  component  of  the  scalar  gradient,  Vr(x.  /  ) , 
will  contribute  to  the  in-i)lane  image-flow  velocity  u.  In  this  civse,  the  equation  for 
the  in-plane  image  flow  can  still  be  written  ;vs. 


Oc  Or 

+  “  ^ 


dc  dc  Or. 


(50) 


where  we  have  jissumed  that  the  image  irradiance  E(x,  f )  is  proportional  to  the  two- 
dimensional  slic(*  of  the  scalar  concentration,  c{x.t),  and  where  the  in-i)lane  im.age- 
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Fig.  12  Contours  of  constant  plane-normal  vorticity.  V  x  u  =  §7  -  •  Solid 

contours  denote  zero-vorticity.  Long  diishes  denote  positive  values  and 
short  (hushes  negative.  Contours  spaced  hy  0.5  s~^ . 

flow  v(docity.  u  =  («.  v}  =  d^/df  (Ecj.  11),  is  the  one  derived  from  the  minimization 
function,  jus  descrilxxl  al)ove. 


Considering  the  transport  of  the  thn>e-dimensional  i.so-scalar  surfac(‘s  (see 
Fig.  15).  we  find  that  the  two-dimensional  u  =  (xi.v)  in-plane  image-flow  v(>locity 
components  arc  related  to  the  threodimensiomil  Uf  =  (uf.  Vf,  W()  fluid-flow  velocity 
and  the  thre(!-diniensional  scalar  gradient  components.  In  particular,  we  have. 


dr. 

= 


OclDx 


Oz  (dc/dx)‘  -f  (dr/Oi/)'^ 


Or 

V  =  V{  +  W{  — 


dc/Oij 


(5i: 


Dz  [dr/Ox)'^  4-  [dejd]))'^ 


As  can  he  seen,  hy  suhstituting  Ecj.  51  in  E(j.  50.  these  relations  recover  the  three- 
dimensional  transj)()rt  e(|uation  for  a  con.serviMl  scalar  field  c(x.t).  in  tin'  c;i.s('  of 
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Fig.  13  Contours  of  constant  divergence,  ^  •  u  =  ^  +  ^.  As  in  Fig.  12.  tin' 
contours  are  spaced  by  0.5  s~'  . 

negligi})le  diffusion,  i.c.. 


Oc  Or  Or  Or  Or.  Or 

— +Uf—  =  ^  +U{—  +V{—  +  W{—  =  0 

Ot  Ox  Ot  Ox  Oy  Oz 


Tlu’se  n^sults  provide  us  with  the  criteria  for  when  tin?  in-plane  image-flow 
velocity  can  be  regarded  ius  a  good  aj)proxiination  to  the  in-phuie  fluid-flow  velocity. 
In  j)articular.  we  will  have. 


and 


U  ~  (if 


)’  ~  I'f  .  if 


—  »  ( 

OcIOz 

\  ^ 

lt'(  ' 

\  (OrlOx)-  +  (OrlOy)-  ^ 

/  Ox 

^  » ( 

(  OrjOz  '' 

\  {drlOx)-  -f  (OrlOy)-  ^ 

1  Oy 

(52b) 
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Fig.  14  Contours  of  constant  second  invariant  of  the  rate-of-strain  tensor  (Eq.  45). 
9-  =  -  5  [  ( i )"  +  ( S )"  -  5  ( W  +  57  )•  ]  ■  Coi.t,n..s  8,.aciKl  by  0.5  . 


Wr*  can  see  that  if  Wf  is  small,  or  if  dcjdz  is  small,  or  both,  by  the  measure  in  Ecj.  52. 
then  the  in-plane  image-flow  velocity  field  can  be  accei)ted  Jis  a  good  representation 
of  the  in-plane  fluid-flow  velocity  field. 

Finally,  since  this  method  estimates  tlie  in-plane  image-flow  velocity  u.  and 
not  the  fluid-flow  velocity  Uf ,  its  application  to  two-dimensional  image  slices  of 
three-dimensional  scalar  field  data  is  identical  to  its  af)plication  to  two-dimensional 
scalar  data. 
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Fig.  15  Two-dimensional  slices  of  a  three-dimensional  scalar  field  c(x,t).  Uf  in¬ 
dicates  the  3-D  fluid-flow  velocity  and  u  is  the  resulting  two-dimensional 
in-plane  image-flow  velocity. 

5.  Voyager  2  images  of  Jupiter 

The  method  is  also  illustrated  on  a  pair  of  the  images  of  the  atmospheric  dy¬ 
namics  of  Jupiter  taken  by  Voyager  2.  These  images  were  taken  from  the  ‘‘Voyager 
Time-Lapse,  Cylindrical-Projection  Jupiter  Mosaics,”  by  Avis  &  Collins  (1983). 
640  X  350  pixel  subimages  of  rotations  349  and  350  were  used  in  the  correlation 
process.  The  subimage  spans  168"  to  97"  lon^tude  and  0"  to  —46"  latitude  (the 
equator  is  at  the  top  of  the  image).  The  subimage  from  rotation  349  is  shown  in 
Fig.  16,  with  an  overlay  of  the  initial  placement  of  the  correlation  volume  neighbor¬ 
hoods.  The  vertical  line  on  the  left  is  a  reference  line  which  is  to  be  deformed  using 
the  mean  zonal  velocities  of  Jupiter  from  Limaye  (1985).  Figure  17  shows  the  same 
region,  one  rotation  later,  with  the  associated  grid  deformed  by  the  nonlinear  cor¬ 
relation  method.  On  the  left  is  the  reference  line  firom  Fig.  16,  carried  the  mean 
zonal  flow.  The  displacement  of  the  centers  of  the  correlation  volumes  is  shown  in 
Fig.  18. 


Fig.  1G  Initial  placement  of  correlation  volumes  ovcrlayed  on  a  sub-image  of  ro¬ 
tation  349  from  the  "Voyager  Time-Lapse,  Cylindrical-Projection  Jiij)iter 
Mosaics.”  Each  square  of  the  grid  denotes  a  correlation  volume.  The  ver¬ 
tical  line  on  the  left  is  a  reference  line  to  be  carried  with  tlu'  mean  zonal 
velocity  of  Jupiter  (see  Fig.  17). 


Fig.  17  Deformation  of  the  correlation  volumes  (see  Fig.  IG).  after  one  rotation. 

The  line  on  the  left  was  deformed  from  the  vertical  line  in  Fig.  IG  using  the 
mean  zonal  velocity  of  Jupiter. 


Fig.  18  Displaceinont  of  grid  control  points  (c<‘nfcrs  of  corrclnt ion  ;ift<  r 

one  rotation.  Tlic  lines  on  tin*  h'ft  denote  the  di>placenieni  mo  the  mean 
zonal  flow  of  .Inpit er. 


6.  Conclusions 


Series  expansions  of  image  displacement,  in  conjinu  lion  with  a  glohal  nouhni’ar 
correlation  method,  can  he  used  to  meiusur<‘  fluid  '••locitie>.  and  vel<>(  ity  gi.idietit' 
from  pairs  of  contiiinons.  convected.  .s<alar  itnages.  It  i>  shi»wn  th.it  iiK  ie.tsmg 
the  order  of  the  expansion  can  improve  the  accuracy  of  the  r«’sult>  I  hi'  propdM  il 
method  does  not  re(inire  di.screte  jiarticles  and  may  also  he  used  in  situations  where 
there  is  a  natural  marker  already  in  the  flow.  i.tj..  sjx'cies  concentration  <an  he 
u.sed  to  mciisun'  velocities  in  compressihli*  flows  T!ie  met  ho. I  i>.  developed  for 
three-dimensional  data  sets  and  demonstrated  on  two-dimensional  images  of  fluid 
flow. 
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